!********************************************************************************** ! This computer software was prepared by Battelle Memorial Institute, hereinafter ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY, ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. ! ! CBMZ module: see module_cbmz.F for references and terms of use !********************************************************************************** module module_cbmz_rodas3_solver contains !----------------------------------------------------------------------- subroutine rodas3_ff_x2( n, t, tnext, hmin, hmax, hstart, & y, abstol, reltol, yposlimit, yneglimit, & yfixed, rconst, & lu_nonzero_v, lu_crow_v, lu_diag_v, lu_icol_v, & info, iok, lunerr, & dydtsubr, yjacsubr, decompsubr, solvesubr ) ! include 'sparse_s.h' ! stiffly accurate rosenbrock 3(2), with ! stiffly accurate embedded formula for error control. ! ! all the arguments aggree with the kpp syntax. ! ! input arguments: ! n = number of dependent variables (i.e., time-varying species) ! [t, tnext] = the integration interval ! hmin, hmax = lower and upper bounds for the selected step-size. ! note that for step = hmin the current computed ! solution is unconditionally accepted by the error ! control mechanism. ! hstart = initial guess step-size ! y = vector of (n) concentrations, contains the ! initial values on input ! abstol, reltol = (n) dimensional vectors of ! componentwise absolute and relative tolerances. ! yposlimit = vector of (n) upper-limit positive concentrations ! yneglimit = vector of (n) upper-limit negative concentrations ! yfixed = vector of (*) concentrations of fixed/non-reacting species ! rconst = vector of (*) concentrations of fixed/non-reacting species ! ! lu_nonzero_v = number of non-zero entries in the "augmented jacobian" ! (i.e., the jacobian with all diagonal elements filled in) ! lu_crow_v = vector of (n) pointers to the crow (?) elements ! lu_diag_v = vector of (n) pointers to the diagonal elements ! of the sparsely organized jacobian. ! jac(lu_diag_v(i)) is the i,i diagonal element ! lu_icol_v = vector of (lu_nonzero_v) pointers to the icol (?) elements ! info(1) = 1 for autonomous system ! = 0 for nonautonomous system ! iok = completion code (output) ! +1 = successful integration ! +2 = successful integration, but some substeps were h=hmin ! and error > abstol,reltol ! -1 = failure -- lu decomposition failed with minimum stepsize h=hmin ! -1001 --> -1999 = failure -- with minimum stepsize h=hmin, ! species i = -(iok+1000) was NaN ! -2001 --> -2999 = failure -- with minimum stepsize h=hmin, ! species i = -(iok+2000) was > yneglimit ! -3001 --> -3999 = failure -- with minimum stepsize h=hmin, ! species i = -(iok+3000) was < yposlimit ! -5001 --> -5999 = failure -- with minimum stepsize h=hmin, ! species i = -(iok+5000) had abs(kn(i)) > ylimit_solvesubr ! before call to solvesubr, where kn is either k1,k2,k3,k4 ! ! dydtsubr = name of routine of derivatives. kpp syntax. ! see the header below. ! yjacsubr = name of routine that computes the jacobian, in ! sparse format. kpp syntax. see the header below. ! decompsubr = name of routine that does sparse lu decomposition ! solvesubr = name of routine that does sparse lu backsolve ! ! output arguments: ! y = the values of concentrations at tend. ! t = equals tend on output. ! info(2) = # of dydtsubr calls. ! info(3) = # of yjacsubr calls. ! info(4) = # of accepted steps. ! info(5) = # of rejected steps. ! info(6) = # of steps that were accepted but had ! (hold.eq.hmin) .and. (err.gt.1) ! which means stepsize=minimum and error>tolerance ! ! adrian sandu, march 1996 ! the center for global and regional environmental research use module_peg_util, only: peg_message, peg_error_fatal implicit none integer nvar_maxd, lu_nonzero_v_maxd parameter (nvar_maxd=99) parameter (lu_nonzero_v_maxd=999) ! common block variables ! subr parameters integer n, info(6), iok, lunerr, lu_nonzero_v integer lu_crow_v(n), lu_diag_v(n), lu_icol_v(lu_nonzero_v) real t, tnext, hmin, hmax, hstart real y(n), abstol(n), reltol(n), yposlimit(n), yneglimit(n) real yfixed(*), rconst(*) external dydtsubr, yjacsubr, decompsubr, solvesubr ! local variables logical isreject, autonom integer nfcn, njac, naccept, nreject, nnocnvg, i, j integer ier real k1(nvar_maxd), k2(nvar_maxd) real k3(nvar_maxd), k4(nvar_maxd) real f1(nvar_maxd), ynew(nvar_maxd) real jac(lu_nonzero_v_maxd) real ghinv, uround real tin, tplus, h, hold, hlowest real err, factor, facmax real dround, c43, tau, x1, x2, ytol real ylimit_solvesubr character*80 errmsg ylimit_solvesubr = 1.0e18 ! check n and lu_nonzero_v if (n .gt. nvar_maxd) then call peg_message( lunerr, '*** rodas3 dimensioning problem' ) write(errmsg,9050) 'n, nvar_maxd = ', n, nvar_maxd call peg_message( lunerr, errmsg ) call peg_error_fatal( lunerr, '*** rodas3 fatal error' ) else if (lu_nonzero_v .gt. lu_nonzero_v_maxd) then call peg_message( lunerr, '*** rodas3 dimensioning problem' ) write(errmsg,9050) 'lu_nonvero_v, lu_nonzero_v_maxd = ', & lu_nonzero_v, lu_nonzero_v_maxd call peg_message( lunerr, errmsg ) call peg_error_fatal( lunerr, '*** rodas3 fatal error' ) end if 9050 format( a, 2(1x,i6) ) ! initialization uround = 1.e-7 dround = sqrt(uround) ! check hmin and hmax hlowest = dround if (info(1) .eq. 1) hlowest = 1.0e-7 if (hmin .lt. hlowest) then call peg_message( lunerr, '*** rodas3 -- hmin is too small' ) write(errmsg,9060) 'hmin and minimum allowed value = ', & hmin, hlowest call peg_message( lunerr, errmsg ) call peg_error_fatal( lunerr, '*** rodas3 fatal error' ) else if (hmin .ge. hmax) then call peg_message( lunerr, '*** rodas3 -- hmin >= hmax' ) write(errmsg,9060) 'hmin, hmax = ', hmin, hmax call peg_message( lunerr, errmsg ) call peg_error_fatal( lunerr, '*** rodas3 fatal error' ) end if 9060 format( a, 1p, 2e14.4 ) ! initialization of counters, etc. autonom = info(1) .eq. 1 !##checkthis## c43 = - 8.e0/3.e0 ! h = 60.e0 ! hmin h = max( hstart, hmin ) tplus = t tin = t isreject = .false. naccept = 0 nreject = 0 nnocnvg = 0 nfcn = 0 njac = 0 ! === starting the time loop === 10 continue tplus = t + h if ( tplus .gt. tnext ) then h = tnext - t tplus = tnext end if call yjacsubr( n, t, y, jac, yfixed, rconst ) njac = njac+1 ghinv = -2.0e0/h do 20 j=1,n jac(lu_diag_v(j)) = jac(lu_diag_v(j)) + ghinv 20 continue call decompsubr( n, jac, ier, lu_crow_v, lu_diag_v, lu_icol_v ) if (ier.ne.0) then if ( h.gt.hmin) then h = 5.0e-1*h go to 10 else ! print *,'ier <> 0, h=',h ! stop iok = -1 goto 200 end if end if call dydtsubr( n, t, y, f1, yfixed, rconst ) ! ====== nonautonomous case =============== if (.not. autonom) then ! tau = sign(dround*max( 1.0e-6, abs(t) ), t) tau = dround*max( 1.0e-6, abs(t) ) call dydtsubr( n, t+tau, y, k2, yfixed, rconst ) nfcn=nfcn+1 do 30 j = 1,n k3(j) = ( k2(j)-f1(j) )/tau 30 continue ! ----- stage 1 (nonautonomous) ----- x1 = 0.5*h do 40 j = 1,n k1(j) = f1(j) + x1*k3(j) if (abs(k1(j)) .gt. ylimit_solvesubr) then iok = -(5000+j) goto 135 end if 40 continue call solvesubr( jac, k1 ) ! ----- stage 2 (nonautonomous) ----- x1 = 4.e0/h x2 = 1.5e0*h do 50 j = 1,n k2(j) = f1(j) - x1*k1(j) + x2*k3(j) if (abs(k2(j)) .gt. ylimit_solvesubr) then iok = -(5000+j) goto 135 end if 50 continue call solvesubr( jac, k2 ) ! ====== autonomous case =============== else ! ----- stage 1 (autonomous) ----- do 60 j = 1,n k1(j) = f1(j) if (abs(k1(j)) .gt. ylimit_solvesubr) then iok = -(5000+j) goto 135 end if 60 continue call solvesubr( jac, k1 ) ! ----- stage 2 (autonomous) ----- x1 = 4.e0/h do 70 j = 1,n k2(j) = f1(j) - x1*k1(j) if (abs(k2(j)) .gt. ylimit_solvesubr) then iok = -(5000+j) goto 135 end if 70 continue call solvesubr( jac, k2 ) end if ! ----- stage 3 ----- do 80 j = 1,n ynew(j) = y(j) - 2.0e0*k1(j) 80 continue call dydtsubr( n, t+h, ynew, f1, yfixed, rconst ) nfcn=nfcn+1 do 90 j = 1,n k3(j) = f1(j) + ( -k1(j) + k2(j) )/h if (abs(k3(j)) .gt. ylimit_solvesubr) then iok = -(5000+j) goto 135 end if 90 continue call solvesubr( jac, k3 ) ! ----- stage 4 ----- do 100 j = 1,n ynew(j) = y(j) - 2.0e0*k1(j) - k3(j) 100 continue call dydtsubr( n, t+h, ynew, f1, yfixed, rconst ) nfcn=nfcn+1 do 110 j = 1,n k4(j) = f1(j) + ( -k1(j) + k2(j) - c43*k3(j) )/h if (abs(k4(j)) .gt. ylimit_solvesubr) then iok = -(5000+j) goto 135 end if 110 continue call solvesubr( jac, k4 ) ! ---- the solution --- do 120 j = 1,n ynew(j) = y(j) - 2.0e0*k1(j) - k3(j) - k4(j) 120 continue ! ====== error estimation ======== err=0.e0 do 130 i=1,n ytol = abstol(i) + reltol(i)*abs(ynew(i)) err = err + ( k4(i)/ytol )**2 130 continue err = max( uround, sqrt( err/n ) ) ! ======= choose the stepsize =============================== factor = 0.9/err**(1.e0/3.e0) if (isreject) then facmax=1.0 else facmax=10.0 end if factor = max( 1.0e-1, min(factor,facmax) ) hold = h h = min( hmax, max(hmin,factor*h) ) ! ======= check for nan, too big, too small ================= ! if any of these conditions occur, either ! exit if hold <= hmin ! reduce step size and retry if hold > hmin iok = 1 do i = 1, n if (y(i) .ne. y(i)) then iok = -(1000+i) goto 135 else if (y(i) .lt. yneglimit(i)) then iok = -(2000+i) goto 135 else if (y(i) .gt. yposlimit(i)) then iok = -(3000+i) goto 135 end if end do 135 if (iok .lt. 0) then if (hold .le. hmin) then goto 200 else isreject = .true. nreject = nreject+1 h = max(hmin, 0.5*hold) if (t.eq.tin) h = max(hmin, 1.0e-1*hold) goto 10 end if end if ! ======= rejected/accepted step ============================ if ( (err.gt.1).and.(hold.gt.hmin) ) then isreject = .true. nreject = nreject+1 if (t.eq.tin) h = max(hmin, 1.0e-1*hold) else isreject = .false. do 140 i=1,n y(i) = ynew(i) 140 continue t = tplus if (err.gt.1) then nnocnvg = nnocnvg+1 else naccept = naccept+1 end if end if ! ======= end of the time loop =============================== if ( t .lt. tnext ) go to 10 iok = 1 if (nnocnvg .gt. 0) iok = 2 ! ======= output information ================================= 200 info(2) = nfcn info(3) = njac info(4) = naccept info(5) = nreject info(6) = nnocnvg return end subroutine rodas3_ff_x2 end module module_cbmz_rodas3_solver