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

   1       SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )
   2
   3       INCLUDE 'KPP_ROOT_params.h'
   4       INCLUDE 'KPP_ROOT_global.h'
   5
   6       INTEGER NSENSIT
   7 C TIN - Start Time
   8       KPP_REAL TIN
   9 C TOUT - End Time
  10       KPP_REAL TOUT
  11 C Y - Concentrations and Sensitivities
  12       KPP_REAL Y(NVAR*(NSENSIT+1))
  13 C ---  Note: Y contains: (1:NVAR) concentrations, followed by
  14 C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
  15 C ---                   etc.,  followed by
  16 C ---                   (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
  17
  18       INTEGER    INFO(5)
  19
  20       EXTERNAL FUNC_CHEM, JAC_CHEM
  21
  22       INFO(1) = Autonomous
  23
  24       CALL ROS2_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,STEPMAX,
  25      +                   STEPMIN,Y,ATOL,RTOL,
  26      +                   Info,FUNC_CHEM,JAC_CHEM)
  27
  28
  29       RETURN
  30       END
  31
  32
  33
  34
  35       SUBROUTINE ROS2_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart,
  36      +                   y,AbsTol,RelTol,
  37      +                   Info,FUNC_CHEM,JAC_CHEM)
  38       IMPLICIT NONE
  39       INCLUDE 'KPP_ROOT_params.h'
  40       INCLUDE 'KPP_ROOT_global.h'
  41       INCLUDE 'KPP_ROOT_sparse.h'
  42 C
  43 C  Ros2 with direct-decoupled calculation of sensitivities
  44 C
  45 C The global variable DDMTYPE distinguishes between:
  46 C      DDMTYPE = 0 : sensitivities w.r.t. initial values
  47 C      DDMTYPE = 1 : sensitivities w.r.t. parameters
  48 C
  49 C  INPUT ARGUMENTS:
  50 C     y = Vector of:   (1:NVAR) concentrations, followed by
  51 C                      (1:NVAR) sensitivities w.r.t. first parameter, followed by
  52 C                       etc.,  followed by
  53 C                      (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
  54 C         (y contains initial values at input, final values at output)
  55 C     [T, Tnext] = the integration interval
  56 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  57 C          Note that for Step = Hmin the current computed
  58 C          solution is unconditionally accepted by the error
  59 C          control mechanism.
  60 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  61 C          componentwise absolute and relative tolerances.
  62 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  63 C          See the header below.
  64 C     JAC_CHEM = name of routine that computes the Jacobian, in
  65 C          sparse format. KPP syntax. See the header below.
  66 C     Info(1) = 1  for  autonomous   system
  67 C             = 0  for nonautonomous system
  68 C
  69 C  OUTPUT ARGUMENTS:
  70 C     y = the values of concentrations at TEND.
  71 C     T = equals TEND on output.
  72 C     Info(2) = # of FUNC_CHEM calls.
  73 C     Info(3) = # of JAC_CHEM calls.
  74 C     Info(4) = # of accepted steps.
  75 C     Info(5) = # of rejected steps.
  76 C
  77 C  Adrian Sandu, December 2001
  78
  79
  80       INTEGER NSENSIT
  81       KPP_REAL y(NVAR*(NSENSIT+1)), ynew(NVAR*(NSENSIT+1))
  82       KPP_REAL K1(NVAR*(NSENSIT+1))
  83       KPP_REAL K2(NVAR*(NSENSIT+1))
  84       KPP_REAL K3(NVAR)
  85       KPP_REAL DFDT(NVAR*(NSENSIT+1))
  86       KPP_REAL DFDP(NVAR*NSENSIT+1), DFDPDT(NVAR*NSENSIT+1)
  87       KPP_REAL DJDP(NVAR*NSENSIT+1)
  88       KPP_REAL F1(NVAR), F2(NVAR)
  89       KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
  90       KPP_REAL DJDT(LU_NONZERO)
  91       KPP_REAL HESS(NHESS)
  92       KPP_REAL Hmin,Hmax,Hnew,Hstart,ghinv,uround
  93       KPP_REAL AbsTol(NVAR), RelTol(NVAR)
  94       KPP_REAL T, Tnext, H, Hold, Tplus, e
  95       KPP_REAL ERR, factor, facmax, dround, elo, tau, gam
  96
  97       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  98       INTEGER    Info(5)
  99       LOGICAL    IsReject,Autonomous,Embed3
 100       EXTERNAL   FUNC_CHEM, JAC_CHEM, HESS_CHEM
 101
 102       LOGICAL    negative
 103       KPP_REAL gamma, m1, m2, alpha, beta, delta, theta, w
 104       KPP_REAL gamma3, d1, d2, d3, beta1, beta2
 105       KPP_REAL c31, c32, c34
 106
 107 c     Initialization of counters, etc.
 108       Autonomous = Info(1) .EQ. 1
 109       Embed3  = Info(2) .EQ. 1
 110       uround = 1.d-15
 111       dround = 1.0d-7 ! DSQRT(uround)
 112       H = DMAX1(1.d-8, Hstart)
 113       Tplus = T
 114       IsReject = .false.
 115       Naccept  = 0
 116       Nreject  = 0
 117       Nfcn     = 0
 118       Njac     = 0
 119       gamma = 1.d0 + 1.d0/DSQRT(2.0d0)
 120       c31 = -1.0D0/gamma**2*(1.0D0-7.0D0*gamma+9.0D0*gamma**2)
 121      &         /(-1.0D0+2.0D0*gamma)
 122       c32 = -1.0D0/gamma**2*(1.0D0-6.0D0*gamma+6.0D0*gamma**2)
 123      &          /(-1.0D0+2.0D0*gamma)/2
 124       gamma3 = 0.5D0 - 2*gamma
 125       d1 = ((-9.0D0*gamma+8.0D0*gamma**2+2.0D0)/gamma**2/
 126      &          (-1.0D0+2*gamma))/6.0D0
 127       d2 = ((-1.0D0+3.0D0*gamma)/gamma**2/
 128      &          (-1.0D0+2.0D0*gamma))/6.0D0
 129       d3 = -1.0D0/(3.0D0*gamma)
 130        m1 = -3.d0/(2.d0*gamma)
 131        m2 = -1.d0/(2.d0*gamma)
 132
 133
 134 C === Starting the time loop ===
 135  10    CONTINUE
 136        Tplus = T + H
 137        IF ( Tplus .gt. Tnext ) THEN
 138           H = Tnext - T
 139           Tplus = Tnext
 140        END IF
 141
 142 C              Initial Function, Jacobian, and Hessian Values
 143        CALL FUNC_CHEM(NVAR, T, y, F1)
 144        CALL JAC_CHEM(NVAR, T, y, JAC)
 145        CALL HESS_CHEM( NVAR, T, y, HESS )
 146        IF (DDMTYPE .EQ. 1) THEN
 147           CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
 148        END IF
 149
 150 C              Estimate the time derivatives in non-autonomous case
 151        IF (.not. Autonomous) THEN
 152          tau = DSIGN(dround*DMAX1( 1.0d0, DABS(T) ), T)
 153          CALL FUNC_CHEM(NVAR, T+tau, y, K2)
 154          nfcn=nfcn+1
 155          CALL JAC_CHEM(NVAR, T+tau, y, AJAC)
 156          njac=njac+1
 157          DO 20 j = 1,NVAR
 158            DFDT(j) = ( K2(j)-F1(j) )/tau
 159  20      CONTINUE
 160          DO 30 j = 1,LU_NONZERO
 161            DJDT(j) = ( AJAC(j)-JAC(j) )/tau
 162  30      CONTINUE
 163          DO 40 i=1,NSENSIT
 164             CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),DFDT(i*NVAR+1))
 165  40      CONTINUE
 166         END IF ! .not. Autonomous
 167
 168        Njac = Njac+1
 169        ghinv = - 1.0d0/(gamma*H)
 170        DO 50 j=1,LU_NONZERO
 171          AJAC(j) = JAC(j)
 172  50    CONTINUE
 173        DO 60 j=1,NVAR
 174          AJAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv
 175  60    CONTINUE
 176        CALL KppDecomp (AJAC, ier)
 177
 178        IF (ier.ne.0) THEN
 179          IF ( H.gt.Hmin) THEN
 180             H = 5.0d-1*H
 181             go to 10
 182          ELSE
 183             print *,'IER <> 0, H=',H
 184             stop
 185          END IF
 186        END IF
 187
 188
 189
 190
 191 C ----- STAGE 1 -----
 192          delta = gamma*H
 193          DO 70 j = 1,NVAR
 194            K1(j) =  F1(j)
 195  70      CONTINUE
 196          IF (.NOT. Autonomous) THEN
 197            DO 80 j = 1,NVAR
 198              K1(j) =  K1(j) + delta*DFDT(j)
 199  80        CONTINUE
 200          END IF
 201          CALL KppSolve (AJAC, K1)
 202 C               --- If  derivative w.r.t. parameters
 203          IF (DDMTYPE .EQ. 1) THEN
 204             CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP)
 205          END IF
 206 C               --- End of derivative w.r.t. parameters
 207          DO 120 i=1,NSENSIT
 208             CALL Jac_SP_Vec (JAC,y(i*NVAR+1),K1(i*NVAR+1))
 209             CALL Hess_Vec ( HESS, K1(1), y(i*NVAR+1), F2 )
 210             DO 90 j=1,NVAR
 211               K1(i*NVAR+j) = K1(i*NVAR+j) + gHinv*F2(j)
 212  90         CONTINUE
 213            IF (.NOT. Autonomous) THEN
 214              DO 100 j = 1,NVAR
 215                K1(i*NVAR+j) =  K1(i*NVAR+j) + delta*DFDT(i*NVAR+j)
 216  100          CONTINUE
 217            END IF
 218 C               --- If  derivative w.r.t. parameters
 219           IF (DDMTYPE .EQ. 1) THEN
 220             DO 110 j = 1,NVAR
 221                K1(i*NVAR+j) = K1(i*NVAR+j) + DFDP((i-1)*NVAR+j)
 222      &                           + DJDP((i-1)*NVAR+j)
 223  110        CONTINUE
 224           END IF
 225 C               --- End of derivative w.r.t. parameters
 226             CALL KppSolve (AJAC, K1(i*NVAR+1))
 227  120     CONTINUE
 228
 229 C ----- STAGE 2   -----
 230        alpha = - 1.d0/gamma
 231        DO 130 j = 1,NVAR*(NSENSIT+1)
 232          ynew(j) = y(j) + alpha*K1(j)
 233  130    CONTINUE
 234        CALL FUNC_CHEM(NVAR, T+H, ynew, F1)
 235        IF (DDMTYPE.EQ.1) THEN
 236          CALL DFUNDPAR(NVAR, NSENSIT, T+H, ynew, DFDP)
 237        END IF
 238        nfcn=nfcn+1
 239        beta1 = 2.d0/(gamma*H)
 240        delta = -gamma*H
 241        DO 140 j = 1,NVAR
 242          K2(j) = F1(j) + beta1*K1(j)
 243  140   CONTINUE
 244        IF (.NOT. Autonomous) THEN
 245           DO 150 j = 1,NVAR
 246                K2(j) =  K2(j) + delta*DFDT(j)
 247  150      CONTINUE
 248        END IF
 249        CALL KppSolve (AJAC, K2)
 250 C               --- If  derivative w.r.t. parameters
 251        IF (DDMTYPE .EQ. 1) THEN
 252             CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP)
 253        END IF
 254 C               --- End of derivative w.r.t. parameters
 255
 256        CALL JAC_CHEM(NVAR, T+H, Ynew, JAC)
 257        njac=njac+1
 258        DO 190 i=1,NSENSIT
 259           CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K2(i*NVAR+1))
 260           CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),F1)
 261           CALL Hess_Vec ( HESS, K2(1), y(i*NVAR+1), F2 )
 262           DO 160 j = 1,NVAR
 263              K2(i*NVAR+j) = K2(i*NVAR+j) + beta1*K1(i*NVAR+j)
 264      &                       + gHinv*F2(j)
 265  160      CONTINUE
 266           IF (.NOT. Autonomous) THEN
 267              DO 170 j = 1,NVAR
 268                K2(i*NVAR+j) =  K2(i*NVAR+j) + delta*DFDT(i*NVAR+j)
 269  170         CONTINUE
 270           END IF
 271 C               --- If  derivative w.r.t. parameters
 272           IF (DDMTYPE .EQ. 1) THEN
 273             DO 180 j = 1,NVAR
 274                K2(i*NVAR+j) = K2(i*NVAR+j)  + DFDP((i-1)*NVAR+j)
 275      &                           + DJDP((i-1)*NVAR+j)
 276  180        CONTINUE
 277           END IF
 278 C               --- End of derivative w.r.t. parameters
 279           CALL KppSolve (AJAC, K2(i*NVAR+1))
 280  190    CONTINUE
 281
 282 C ----- STAGE 3  for error control only -----
 283        IF (Embed3) THEN
 284        beta1 = -c31/H
 285        beta2 = -c32/H
 286        delta = gamma3*H
 287        DO 195 j = 1,NVAR
 288          K3(j) = F1(j) + beta1*K1(j) + beta2*K2(j)
 289  195   CONTINUE
 290        IF (.NOT. Autonomous) THEN
 291          DO 196 j = 1,NVAR
 292            K3(j) = K3(j) + delta*DFDT(j)
 293  196     CONTINUE
 294        END IF
 295        CALL KppSolve (AJAC, K3)
 296        END IF
 297
 298 C ---- The Solution ---
 299        DO 200 j = 1,NVAR*(NSENSIT+1)
 300          ynew(j) = y(j) + m1*K1(j) + m2*K2(j)
 301  200   CONTINUE
 302
 303
 304 C ====== Error estimation for concentrations only; this can be easily adapted to
 305 C                                             estimate the sensitivity error too ========
 306
 307         ERR=0.d0
 308         DO 210 i=1,NVAR
 309            w = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(ynew(i)))
 310            IF (Embed3) THEN
 311              e = d1*K1(i) + d2*K2(i) + d3*K3(i)
 312            ELSE
 313              e = (1.d0/(2.d0*gamma))*(K1(i)+K2(i))
 314            END IF
 315            ERR = ERR + ( e/w )**2
 316  210    CONTINUE
 317         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 318
 319 C ======= Choose the stepsize ===============================
 320
 321         IF (Embed3) THEN
 322            elo    = 3.0D0 ! estimator local order
 323         ELSE
 324            elo    = 2.0D0
 325         END IF
 326         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 327         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 328
 329 C ======= Rejected/Accepted Step ============================
 330
 331         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 332           IsReject = .true.
 333           H = DMIN1(H/10,Hnew)
 334           Nreject  = Nreject+1
 335         ELSE
 336           DO 300 i=1,NVAR*(NSENSIT+1)
 337              y(i)  = ynew(i)
 338  300      CONTINUE
 339           T = Tplus
 340           IF (.NOT.IsReject) THEN
 341               H = Hnew   ! Do not increase stepsize if previous step was rejected
 342           END IF
 343           IsReject = .false.
 344           Naccept = Naccept+1
 345         END IF
 346
 347 C ======= End of the time loop ===============================
 348       IF ( T .lt. Tnext ) GO TO 10
 349
 350
 351
 352 C ======= Output Information =================================
 353       Info(2) = Nfcn
 354       Info(3) = Njac
 355       Info(4) = Naccept
 356       Info(5) = Nreject
 357
 358       RETURN
 359       END
 360
 361
 362
 363       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 364       INCLUDE 'KPP_ROOT_params.h'
 365       INCLUDE 'KPP_ROOT_global.h'
 366       KPP_REAL   T, Told
 367       KPP_REAL   Y(NVAR), P(NVAR)
 368       Told = TIME
 369       TIME = T
 370       CALL Update_SUN()
 371       CALL Update_RCONST()
 372       CALL Fun( Y,  FIX, RCONST, P )
 373       TIME = Told
 374       RETURN
 375       END
 376
 377
 378       SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
 379 C ---  Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
 380       INCLUDE 'KPP_ROOT_params.h'
 381       INCLUDE 'KPP_ROOT_global.h'
 382 C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
 383       INTEGER NCOEFF, JCOEFF(NREACT)
 384       COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
 385
 386       KPP_REAL   T, Told
 387       KPP_REAL   Y(NVAR), P(NVAR*NSENSIT)
 388       Told = TIME
 389       TIME = T
 390       CALL Update_SUN()
 391       CALL Update_RCONST()
 392 C
 393       IF (DDMTYPE .EQ. 0) THEN
 394 C ---  Note: the values below are for sensitivities w.r.t. initial values;
 395 C ---       they may have to be changed for other applications
 396         DO j=1,NSENSIT
 397           DO i=1,NVAR
 398             P(i+NVAR*(j-1)) = 0.0D0
 399           END DO
 400         END DO
 401       ELSE
 402 C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
 403 C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
 404 C ---       w.r.t. which one differentiates
 405         CALL dFun_dRcoeff( Y,  FIX, NCOEFF, JCOEFF, P )
 406       END IF
 407       TIME = Told
 408       RETURN
 409       END
 410
 411       SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
 412 C ---  Computes the partial derivatives of JAC w.r.t. parameters times user vector U
 413       INCLUDE 'KPP_ROOT_params.h'
 414       INCLUDE 'KPP_ROOT_global.h'
 415 C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
 416       INTEGER NCOEFF, JCOEFF(NREACT)
 417       COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
 418
 419       KPP_REAL   T, Told
 420       KPP_REAL   Y(NVAR), U(NVAR)
 421       KPP_REAL   P(NVAR*NSENSIT)
 422       Told = TIME
 423       TIME = T
 424       CALL Update_SUN()
 425       CALL Update_RCONST()
 426 C
 427       IF (DDMTYPE .EQ. 0) THEN
 428 C ---  Note: the values below are for sensitivities w.r.t. initial values;
 429 C ---       they may have to be changed for other applications
 430         DO j=1,NSENSIT
 431           DO i=1,NVAR
 432             P(i+NVAR*(j-1)) = 0.0D0
 433           END DO
 434         END DO
 435       ELSE
 436 C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
 437 C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
 438 C ---       w.r.t. which one differentiates
 439         CALL dJac_dRcoeff( Y,  FIX, U, NCOEFF, JCOEFF, P )
 440       END IF
 441       TIME = Told
 442       RETURN
 443       END
 444
 445       SUBROUTINE JAC_CHEM(N, T, Y, J)
 446       INCLUDE 'KPP_ROOT_params.h'
 447       INCLUDE 'KPP_ROOT_global.h'
 448       INTEGER N
 449       KPP_REAL   Told, T
 450       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 451       Told = TIME
 452       TIME = T
 453       CALL Update_SUN()
 454       CALL Update_RCONST()
 455       CALL Jac_SP( Y, FIX, RCONST, J )
 456       TIME = Told
 457       RETURN
 458       END
 459
 460
 461       SUBROUTINE HESS_CHEM(N, T, Y, HESS)
 462       INCLUDE 'KPP_ROOT_params.h'
 463       INCLUDE 'KPP_ROOT_global.h'
 464       INTEGER N
 465       KPP_REAL   Told, T
 466       KPP_REAL   Y(NVAR), HESS(NHESS)
 467       Told = TIME
 468       TIME = T
 469       CALL Update_SUN()
 470       CALL Update_RCONST()
 471       CALL Hessian( Y, FIX, RCONST, HESS )
 472       TIME = Told
 473       RETURN
 474       END
 475
 476
 477
 478
 479
 480