chem/KPP/kpp/kpp-2.1/int/oldies/ros1_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 TOUT - End Time
  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,  HESS_CHEM
  21
  22       INFO(1) = Autonomous
  23
  24       CALL ROS1_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,Y,
  25      +                   Info,FUNC_CHEM,JAC_CHEM,HESS_CHEM)
  26
  27
  28       RETURN
  29       END
  30
  31
  32
  33
  34       SUBROUTINE ROS1_DDM(N,NSENSIT,T,Tnext,Hstart,
  35      +                   y,Info,FUNC_CHEM,JAC_CHEM,HESS_CHEM)
  36       IMPLICIT NONE
  37       INCLUDE 'KPP_ROOT_params.h'
  38       INCLUDE 'KPP_ROOT_sparse.h'
  39       INCLUDE 'KPP_ROOT_global.h'
  40 C
  41 C  Linearly Implicit Euler with direct-decoupled calculation of sensitivities
  42 C  A method of theoretical interest but of no practical value
  43 C
  44 C The global variable DDMTYPE distinguishes between:
  45 C      DDMTYPE = 0 : sensitivities w.r.t. initial values
  46 C      DDMTYPE = 1 : sensitivities w.r.t. parameters
  47 C
  48 C  INPUT ARGUMENTS:
  49 C     y = Vector of:   (1:NVAR) concentrations, followed by
  50 C                      (1:NVAR) sensitivities w.r.t. first parameter, followed by
  51 C                       etc.,  followed by
  52 C                      (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
  53 C         (y contains initial values at input, final values at output)
  54 C     [T, Tnext] = the integration interval
  55 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  56 C          Note that for Step = Hmin the current computed
  57 C          solution is unconditionally accepted by the error
  58 C          control mechanism.
  59 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  60 C          componentwise absolute and relative tolerances.
  61 C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
  62 C          See the header below.
  63 C     JAC_CHEM = name of routine that computes the Jacobian, in
  64 C          sparse format. KPP syntax. See the header below.
  65 C     Info(1) = 1  for  Autonomous   system
  66 C             = 0  for nonAutonomous system
  67 C
  68 C  OUTPUT ARGUMENTS:
  69 C     y = the values of concentrations at Tend.
  70 C     T = equals TENDon output.
  71 C     Info(2) = # of FUNC_CHEM CALLs.
  72 C     Info(3) = # of JAC_CHEM CALLs.
  73 C     Info(4) = # of accepted steps.
  74 C     Info(5) = # of rejected steps.
  75 C     Hstart = The last accepted stepsize
  76 C
  77 C    Adrian Sandu, December 2001
  78 C
  79       INTEGER   NSENSIT
  80       KPP_REAL Fv(NVAR*(NSENSIT+1)), Hv(NVAR)
  81       KPP_REAL DFDP(NVAR*NSENSIT)
  82       KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
  83       KPP_REAL HESS(NHESS)
  84       KPP_REAL DJDP(NVAR*NSENSIT)
  85       KPP_REAL H, Hstart
  86       KPP_REAL y(NVAR*(NSENSIT+1))
  87       KPP_REAL T, Tnext, Tplus
  88       KPP_REAL elo,ghinv,uround
  89
  90       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j,ier
  91       INTEGER    Info(5)
  92       LOGICAL    IsReject, Autonomous
  93       EXTERNAL   FUNC_CHEM, JAC_CHEM, HESS_CHEM
  94
  95
  96       H     = Hstart
  97       Tplus = T
  98       Nfcn  = 0
  99       Njac  = 0
 100
 101 C === Starting the time loop ===
 102  10    CONTINUE
 103
 104        Tplus = T + H
 105        IF ( Tplus .gt. Tnext ) THEN
 106           H = Tnext - T
 107           Tplus = Tnext
 108        END IF
 109
 110 C              Initial Function and Jacobian values
 111        CALL FUNC_CHEM(NVAR, T, y, Fv)
 112        Nfcn = Nfcn+1
 113        CALL JAC_CHEM(NVAR, T, y, JAC)
 114        Njac = Njac+1
 115        CALL HESS_CHEM( NVAR, T, y, HESS )
 116        IF (DDMTYPE .EQ. 1) THEN
 117           CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
 118        END IF
 119
 120 C              Form the Prediction matrix and compute its LU factorization
 121        DO 40 j=1,LU_NONZERO
 122          AJAC(j) = -JAC(j)
 123  40    CONTINUE
 124        DO 50 j=1,NVAR
 125          AJAC(LU_DIAG(j)) = AJAC(LU_DIAG(j)) + 1.0d0/H
 126  50    CONTINUE
 127        CALL KppDecomp (AJAC, ier)
 128 C
 129        IF (ier.ne.0) THEN
 130           PRINT *,'ROS1: Singular factorization at T=',T,'; H=',H
 131           STOP
 132        END IF
 133
 134 C ------------ STAGE 1-------------------------
 135        CALL KppSolve (AJAC, Fv)
 136 C               --- If  derivative w.r.t. parameters
 137          IF (DDMTYPE .EQ. 1) THEN
 138             CALL DJACDPAR(NVAR, NSENSIT, T, y, Fv(1), DJDP)
 139          END IF
 140 C               --- End of derivative w.r.t. parameters
 141
 142        DO 100 i=1,NSENSIT
 143           CALL Jac_SP_Vec (JAC, y(i*NVAR+1), Fv(i*NVAR+1))
 144           CALL Hess_Vec ( HESS, y(i*NVAR+1), Fv(1), Hv )
 145           IF (DDMTYPE .EQ. 0) THEN
 146             DO 80 j=1,NVAR
 147               Fv(i*NVAR+j) = Fv(i*NVAR+j) + Hv(j)
 148  80         CONTINUE
 149           ELSE
 150             DO 90 j=1,NVAR
 151               Fv(i*NVAR+j) = Fv(i*NVAR+j) + Hv(j)
 152      &                + DFDP(i*NVAR+j)+ DJDP((i-1)*NVAR+j)
 153  90         CONTINUE
 154           END IF
 155           CALL KppSolve (AJAC, Fv(i*NVAR+1))
 156  100   CONTINUE
 157
 158 C ---- The Solution ---
 159        DO 160 j = 1,NVAR*(NSENSIT+1)
 160          y(j) = y(j) + Fv(j)
 161  160   CONTINUE
 162        T = T + H
 163
 164 C ======= End of the time loop ===============================
 165       IF ( T .lt. Tnext ) THEN
 166           GO TO 10
 167       END IF
 168
 169 C ======= Output Information =================================
 170       Info(2) = Nfcn
 171       Info(3) = Njac
 172       Info(4) = Naccept
 173       Info(5) = Nreject
 174       Hstart  = H
 175
 176       RETURN
 177       END
 178
 179
 180
 181       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 182       INCLUDE 'KPP_ROOT_params.h'
 183       INCLUDE 'KPP_ROOT_global.h'
 184       INTEGER N
 185       KPP_REAL   T, Told
 186       KPP_REAL   Y(NVAR), P(NVAR)
 187       Told = TIME
 188       TIME = T
 189       CALL Update_SUN()
 190       CALL Update_RCONST()
 191       CALL Fun( Y,  FIX, RCONST, P )
 192       TIME = Told
 193       RETURN
 194       END
 195
 196
 197       SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
 198 C ---  Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
 199       INCLUDE 'KPP_ROOT_params.h'
 200       INCLUDE 'KPP_ROOT_global.h'
 201 C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
 202       INTEGER NCOEFF, JCOEFF(NREACT)
 203       COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
 204
 205       INTEGER N
 206       KPP_REAL   T, Told
 207       KPP_REAL   Y(NVAR), P(NVAR*NSENSIT)
 208       Told = TIME
 209       TIME = T
 210       CALL Update_SUN()
 211       CALL Update_RCONST()
 212 C
 213       IF (DDMTYPE .EQ. 0) THEN
 214 C ---  Note: the values below are for sensitivities w.r.t. initial values;
 215 C ---       they may have to be changed for other applications
 216         DO j=1,NSENSIT
 217           DO i=1,NVAR
 218             P(i+NVAR*(j-1)) = 0.0D0
 219           END DO
 220         END DO
 221       ELSE
 222 C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
 223 C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
 224 C ---       w.r.t. which one differentiates
 225         CALL dFun_dRcoeff( Y,  FIX, NCOEFF, JCOEFF, P )
 226       END IF
 227       TIME = Told
 228       RETURN
 229       END
 230
 231       SUBROUTINE JAC_CHEM(N, T, Y, J)
 232       INCLUDE 'KPP_ROOT_params.h'
 233       INCLUDE 'KPP_ROOT_global.h'
 234       INTEGER N
 235       KPP_REAL   Told, T
 236       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 237       Told = TIME
 238       TIME = T
 239       CALL Update_SUN()
 240       CALL Update_RCONST()
 241       CALL Jac_SP( Y,  FIX, RCONST, J )
 242       TIME = Told
 243       RETURN
 244       END
 245
 246
 247       SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
 248 C ---  Computes the partial derivatives of JAC w.r.t. parameters times user vector U
 249       INCLUDE 'KPP_ROOT_params.h'
 250       INCLUDE 'KPP_ROOT_global.h'
 251 C ---  NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
 252       INTEGER NCOEFF, JCOEFF(NREACT)
 253       COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
 254
 255       INTEGER N
 256       KPP_REAL   T, Told
 257       KPP_REAL   Y(NVAR), U(NVAR)
 258       KPP_REAL   P(NVAR*NSENSIT)
 259       Told = TIME
 260       TIME = T
 261       CALL Update_SUN()
 262       CALL Update_RCONST()
 263 C
 264       IF (DDMTYPE .EQ. 0) THEN
 265 C ---  Note: the values below are for sensitivities w.r.t. initial values;
 266 C ---       they may have to be changed for other applications
 267         DO j=1,NSENSIT
 268           DO i=1,NVAR
 269             P(i+NVAR*(j-1)) = 0.0D0
 270           END DO
 271         END DO
 272       ELSE
 273 C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
 274 C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
 275 C ---       w.r.t. which one differentiates
 276         CALL dJac_dRcoeff( Y,  FIX, U, NCOEFF, JCOEFF, P )
 277       END IF
 278       TIME = Told
 279       RETURN
 280       END
 281
 282
 283       SUBROUTINE HESS_CHEM(N, T, Y, HESS)
 284       INCLUDE 'KPP_ROOT_params.h'
 285       INCLUDE 'KPP_ROOT_global.h'
 286       INTEGER N
 287       KPP_REAL   Told, T
 288       KPP_REAL   Y(NVAR), HESS(NHESS)
 289       Told = TIME
 290       TIME = T
 291       CALL Update_SUN()
 292       CALL Update_RCONST()
 293       CALL Hessian( Y, FIX, RCONST, HESS )
 294       TIME = Told
 295       RETURN
 296       END
 297
 298
 299
 300