chem/module_cbmz_rodas3_solver.F

   1 !**********************************************************************************
   2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
   3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
   4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
   5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
   6 !
   7 ! CBMZ module: see module_cbmz.F for references and terms of use
   8 !**********************************************************************************
   9       module module_cbmz_rodas3_solver
  10
  11       contains
  12
  13
  14 !-----------------------------------------------------------------------
  15
  16       subroutine rodas3_ff_x2( n, t, tnext, hmin, hmax, hstart,   &
  17             y, abstol, reltol, yposlimit, yneglimit,   &
  18             yfixed, rconst,   &
  19             lu_nonzero_v, lu_crow_v, lu_diag_v, lu_icol_v,   &
  20             info, iok, lunerr,   &
  21             dydtsubr, yjacsubr, decompsubr, solvesubr )
  22 !     include 'sparse_s.h'
  23
  24 !       stiffly accurate rosenbrock 3(2), with
  25 !     stiffly accurate embedded formula for error control.
  26 !
  27 !     all the arguments aggree with the kpp syntax.
  28 !
  29 !  input arguments:
  30 !     n = number of dependent variables (i.e., time-varying species)
  31 !     [t, tnext] = the integration interval
  32 !     hmin, hmax = lower and upper bounds for the selected step-size.
  33 !          note that for step = hmin the current computed
  34 !          solution is unconditionally accepted by the error
  35 !          control mechanism.
  36 !     hstart = initial guess step-size
  37 !     y = vector of (n) concentrations, contains the
  38 !         initial values on input
  39 !     abstol, reltol = (n) dimensional vectors of
  40 !          componentwise absolute and relative tolerances.
  41 !     yposlimit = vector of (n) upper-limit positive concentrations
  42 !     yneglimit = vector of (n) upper-limit negative concentrations
  43 !     yfixed = vector of (*) concentrations of fixed/non-reacting species
  44 !     rconst = vector of (*) concentrations of fixed/non-reacting species
  45 !
  46 !     lu_nonzero_v = number of non-zero entries in the "augmented jacobian"
  47 !          (i.e., the jacobian with all diagonal elements filled in)
  48 !     lu_crow_v = vector of (n) pointers to the crow (?) elements
  49 !     lu_diag_v = vector of (n) pointers to the diagonal elements
  50 !          of the sparsely organized jacobian.
  51 !          jac(lu_diag_v(i)) is the i,i diagonal element
  52 !     lu_icol_v = vector of (lu_nonzero_v) pointers to the icol (?) elements
  53 !     info(1) = 1  for  autonomous   system
  54 !             = 0  for nonautonomous system
  55 !     iok = completion code (output)
  56 !       +1 = successful integration
  57 !       +2 = successful integration, but some substeps were h=hmin
  58 !               and error > abstol,reltol
  59 !       -1 = failure -- lu decomposition failed with minimum stepsize h=hmin
  60 !       -1001 --> -1999 = failure -- with minimum stepsize h=hmin,
  61 !               species i = -(iok+1000) was NaN
  62 !       -2001 --> -2999 = failure -- with minimum stepsize h=hmin,
  63 !               species i = -(iok+2000) was > yneglimit
  64 !       -3001 --> -3999 = failure -- with minimum stepsize h=hmin,
  65 !               species i = -(iok+3000) was < yposlimit
  66 !       -5001 --> -5999 = failure -- with minimum stepsize h=hmin,
  67 !               species i = -(iok+5000) had abs(kn(i)) > ylimit_solvesubr
  68 !               before call to solvesubr, where kn is either k1,k2,k3,k4
  69 !
  70 !     dydtsubr = name of routine of derivatives. kpp syntax.
  71 !          see the header below.
  72 !     yjacsubr = name of routine that computes the jacobian, in
  73 !          sparse format. kpp syntax. see the header below.
  74 !     decompsubr = name of routine that does sparse lu decomposition
  75 !     solvesubr = name of routine that does sparse lu backsolve
  76 !
  77 !  output arguments:
  78 !     y = the values of concentrations at tend.
  79 !     t = equals tend on output.
  80 !     info(2) = # of dydtsubr calls.
  81 !     info(3) = # of yjacsubr calls.
  82 !     info(4) = # of accepted steps.
  83 !     info(5) = # of rejected steps.
  84 !     info(6) = # of steps that were accepted but had
  85 !               (hold.eq.hmin) .and. (err.gt.1)
  86 !               which means stepsize=minimum and error>tolerance
  87 !
  88 !     adrian sandu, march 1996
  89 !     the center for global and regional environmental research
  90
  91       use module_peg_util, only:  peg_message, peg_error_fatal
  92
  93       implicit none
  94
  95       integer nvar_maxd, lu_nonzero_v_maxd
  96       parameter (nvar_maxd=99)
  97       parameter (lu_nonzero_v_maxd=999)
  98
  99 ! common block variables
 100
 101 ! subr parameters
 102       integer    n, info(6), iok, lunerr, lu_nonzero_v
 103       integer    lu_crow_v(n), lu_diag_v(n), lu_icol_v(lu_nonzero_v)
 104       real       t, tnext, hmin, hmax, hstart
 105       real       y(n), abstol(n), reltol(n), yposlimit(n), yneglimit(n)
 106       real       yfixed(*), rconst(*)
 107       external   dydtsubr, yjacsubr, decompsubr, solvesubr
 108
 109 ! local variables
 110       logical    isreject, autonom
 111       integer    nfcn, njac, naccept, nreject, nnocnvg, i, j
 112       integer    ier
 113
 114       real       k1(nvar_maxd), k2(nvar_maxd)
 115       real       k3(nvar_maxd), k4(nvar_maxd)
 116       real       f1(nvar_maxd), ynew(nvar_maxd)
 117       real       jac(lu_nonzero_v_maxd)
 118       real       ghinv, uround
 119       real       tin, tplus, h, hold, hlowest
 120       real       err, factor, facmax
 121       real       dround, c43, tau, x1, x2, ytol
 122       real       ylimit_solvesubr
 123
 124       character*80 errmsg
 125
 126       ylimit_solvesubr = 1.0e18
 127
 128 !     check n and lu_nonzero_v
 129       if (n .gt. nvar_maxd) then
 130           call peg_message( lunerr, '*** rodas3 dimensioning problem' )
 131           write(errmsg,9050) 'n, nvar_maxd = ', n, nvar_maxd
 132           call peg_message( lunerr, errmsg )
 133           call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
 134       else if (lu_nonzero_v .gt. lu_nonzero_v_maxd) then
 135           call peg_message( lunerr, '*** rodas3 dimensioning problem' )
 136           write(errmsg,9050) 'lu_nonvero_v, lu_nonzero_v_maxd = ',   &
 137                 lu_nonzero_v, lu_nonzero_v_maxd
 138           call peg_message( lunerr, errmsg )
 139           call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
 140       end if
 141 9050  format( a, 2(1x,i6) )
 142
 143 !     initialization
 144       uround = 1.e-7
 145       dround = sqrt(uround)
 146
 147 !     check hmin and hmax
 148       hlowest = dround
 149       if (info(1) .eq. 1) hlowest = 1.0e-7
 150       if (hmin .lt. hlowest) then
 151           call peg_message( lunerr, '*** rodas3 -- hmin is too small' )
 152           write(errmsg,9060) 'hmin and minimum allowed value = ',   &
 153                 hmin, hlowest
 154           call peg_message( lunerr, errmsg )
 155           call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
 156       else if (hmin .ge. hmax) then
 157           call peg_message( lunerr, '*** rodas3 -- hmin >= hmax' )
 158           write(errmsg,9060) 'hmin, hmax = ', hmin, hmax
 159           call peg_message( lunerr, errmsg )
 160           call peg_error_fatal( lunerr, '*** rodas3 fatal error' )
 161       end if
 162 9060  format( a, 1p, 2e14.4 )
 163
 164 !     initialization of counters, etc.
 165       autonom = info(1) .eq. 1
 166 !##checkthis##
 167       c43 = - 8.e0/3.e0
 168 !     h = 60.e0 ! hmin
 169       h = max( hstart, hmin )
 170       tplus = t
 171       tin = t
 172       isreject = .false.
 173       naccept  = 0
 174       nreject  = 0
 175       nnocnvg  = 0
 176       nfcn     = 0
 177       njac     = 0
 178
 179
 180 ! === starting the time loop ===
 181  10    continue
 182        tplus = t + h
 183        if ( tplus .gt. tnext ) then
 184           h = tnext - t
 185           tplus = tnext
 186        end if
 187
 188        call yjacsubr( n, t, y, jac, yfixed, rconst )
 189        njac = njac+1
 190        ghinv = -2.0e0/h
 191        do 20 j=1,n
 192          jac(lu_diag_v(j)) = jac(lu_diag_v(j)) + ghinv
 193  20    continue
 194        call decompsubr( n, jac, ier, lu_crow_v, lu_diag_v, lu_icol_v )
 195
 196        if (ier.ne.0) then
 197          if ( h.gt.hmin) then
 198             h = 5.0e-1*h
 199             go to 10
 200          else
 201 !           print *,'ier <> 0, h=',h
 202 !           stop
 203             iok = -1
 204             goto 200
 205          end if
 206        end if
 207
 208        call dydtsubr( n, t, y, f1, yfixed, rconst )
 209
 210 ! ====== nonautonomous case ===============
 211        if (.not. autonom) then
 212 !        tau = sign(dround*max( 1.0e-6, abs(t) ), t)
 213          tau = dround*max( 1.0e-6, abs(t) )
 214          call dydtsubr( n, t+tau, y, k2, yfixed, rconst )
 215          nfcn=nfcn+1
 216          do 30 j = 1,n
 217            k3(j) = ( k2(j)-f1(j) )/tau
 218  30      continue
 219
 220 ! ----- stage 1 (nonautonomous) -----
 221          x1 = 0.5*h
 222          do 40 j = 1,n
 223            k1(j) =  f1(j) + x1*k3(j)
 224            if (abs(k1(j)) .gt. ylimit_solvesubr) then
 225                iok = -(5000+j)
 226                goto 135
 227            end if
 228  40      continue
 229          call solvesubr( jac, k1 )
 230
 231 ! ----- stage 2 (nonautonomous) -----
 232          x1 = 4.e0/h
 233          x2 = 1.5e0*h
 234          do 50 j = 1,n
 235            k2(j) = f1(j) - x1*k1(j) + x2*k3(j)
 236            if (abs(k2(j)) .gt. ylimit_solvesubr) then
 237                iok = -(5000+j)
 238                goto 135
 239            end if
 240  50      continue
 241          call solvesubr( jac, k2 )
 242
 243 ! ====== autonomous case ===============
 244        else
 245 ! ----- stage 1 (autonomous) -----
 246          do 60 j = 1,n
 247            k1(j) =  f1(j)
 248            if (abs(k1(j)) .gt. ylimit_solvesubr) then
 249                iok = -(5000+j)
 250                goto 135
 251            end if
 252  60      continue
 253          call solvesubr( jac, k1 )
 254
 255 ! ----- stage 2 (autonomous) -----
 256          x1 = 4.e0/h
 257          do 70 j = 1,n
 258            k2(j) = f1(j) - x1*k1(j)
 259            if (abs(k2(j)) .gt. ylimit_solvesubr) then
 260                iok = -(5000+j)
 261                goto 135
 262            end if
 263  70      continue
 264          call solvesubr( jac, k2 )
 265        end if
 266
 267 ! ----- stage 3 -----
 268        do 80 j = 1,n
 269          ynew(j) = y(j) - 2.0e0*k1(j)
 270  80    continue
 271        call dydtsubr( n, t+h, ynew, f1, yfixed, rconst )
 272        nfcn=nfcn+1
 273        do 90 j = 1,n
 274          k3(j) = f1(j) + ( -k1(j) + k2(j) )/h
 275          if (abs(k3(j)) .gt. ylimit_solvesubr) then
 276              iok = -(5000+j)
 277              goto 135
 278          end if
 279  90    continue
 280        call solvesubr( jac, k3 )
 281
 282 ! ----- stage 4 -----
 283        do 100 j = 1,n
 284          ynew(j) = y(j) - 2.0e0*k1(j) - k3(j)
 285  100   continue
 286        call dydtsubr( n, t+h, ynew, f1, yfixed, rconst )
 287        nfcn=nfcn+1
 288        do 110 j = 1,n
 289          k4(j) = f1(j) + ( -k1(j) + k2(j) - c43*k3(j)  )/h
 290          if (abs(k4(j)) .gt. ylimit_solvesubr) then
 291              iok = -(5000+j)
 292              goto 135
 293          end if
 294  110   continue
 295        call solvesubr( jac, k4 )
 296
 297 ! ---- the solution ---
 298
 299        do 120 j = 1,n
 300          ynew(j) = y(j) - 2.0e0*k1(j) - k3(j) - k4(j)
 301  120   continue
 302
 303
 304 ! ====== error estimation ========
 305
 306         err=0.e0
 307         do 130 i=1,n
 308            ytol = abstol(i) + reltol(i)*abs(ynew(i))
 309            err = err + ( k4(i)/ytol )**2
 310  130    continue
 311         err = max( uround, sqrt( err/n ) )
 312
 313 ! ======= choose the stepsize ===============================
 314
 315         factor = 0.9/err**(1.e0/3.e0)
 316         if (isreject) then
 317             facmax=1.0
 318         else
 319             facmax=10.0
 320         end if
 321         factor = max( 1.0e-1, min(factor,facmax) )
 322         hold = h
 323         h = min( hmax, max(hmin,factor*h) )
 324
 325 ! ======= check for nan, too big, too small =================
 326 !   if any of these conditions occur, either
 327 !       exit if hold <= hmin
 328 !       reduce step size and retry if hold > hmin
 329
 330         iok = 1
 331         do i = 1, n
 332             if (y(i) .ne. y(i)) then
 333                 iok = -(1000+i)
 334                 goto 135
 335             else if (y(i) .lt. yneglimit(i)) then
 336                 iok = -(2000+i)
 337                 goto 135
 338             else if (y(i) .gt. yposlimit(i)) then
 339                 iok = -(3000+i)
 340                 goto 135
 341             end if
 342         end do
 343 135     if (iok .lt. 0) then
 344             if (hold .le. hmin) then
 345                 goto 200
 346             else
 347                 isreject = .true.
 348                 nreject  = nreject+1
 349                 h = max(hmin, 0.5*hold)
 350                 if (t.eq.tin) h = max(hmin, 1.0e-1*hold)
 351                 goto 10
 352             end if
 353         end if
 354
 355
 356 ! ======= rejected/accepted step ============================
 357
 358         if ( (err.gt.1).and.(hold.gt.hmin) ) then
 359           isreject = .true.
 360           nreject  = nreject+1
 361           if (t.eq.tin) h = max(hmin, 1.0e-1*hold)
 362         else
 363           isreject = .false.
 364           do 140 i=1,n
 365              y(i)  = ynew(i)
 366  140      continue
 367           t = tplus
 368           if (err.gt.1) then
 369             nnocnvg = nnocnvg+1
 370           else
 371             naccept = naccept+1
 372           end if
 373         end if
 374
 375 ! ======= end of the time loop ===============================
 376       if ( t .lt. tnext ) go to 10
 377
 378       iok = 1
 379       if (nnocnvg .gt. 0) iok = 2
 380
 381
 382 ! ======= output information =================================
 383 200   info(2) = nfcn
 384       info(3) = njac
 385       info(4) = naccept
 386       info(5) = nreject
 387       info(6) = nnocnvg
 388
 389       return
 390       end subroutine rodas3_ff_x2
 391
 392
 393       end module module_cbmz_rodas3_solver