[2696] | 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 | |
---|