!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! SDIRK-TLM - Tangent Linear Model of Singly-Diagonally-Implicit RK ! ! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 ! ! * Sdirk 3a: L-stable, 3 stages, order 2, adj-invariant ! ! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! ! (C) Adrian Sandu, July 2005 ! ! Virginia Polytechnic Institute and State University ! ! Contact: sandu@cs.vt.edu ! ! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! ! This implementation is part of KPP - the Kinetic PreProcessor ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! MODULE KPP_ROOT_Integrator USE KPP_ROOT_Precision USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG USE KPP_ROOT_Jacobian, ONLY: Jac_SP_Vec USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, & KppSolve, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD IMPLICIT NONE PUBLIC SAVE !~~~> Statistics on the work performed by the SDIRK method INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, & Nrej=5, Ndec=6, Nsol=7, Nsng=8, & Ntexit=1, Nhexit=2, Nhnew=3 CONTAINS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE INTEGRATE_TLM( NTLM, Y, Y_tlm, TIN, TOUT, ATOL_tlm,RTOL_tlm, & ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE KPP_ROOT_Parameters, ONLY: NVAR,ind_O3 USE KPP_ROOT_Global, ONLY: ATOL,RTOL,VAR IMPLICIT NONE !~~~> Y - Concentrations KPP_REAL :: Y(NVAR) !~~~> NTLM - No. of sensitivity coefficients INTEGER NTLM !~~~> Y_tlm - Sensitivities of concentrations ! Note: Y_tlm (1:NVAR,j) contains sensitivities of ! Y(1:NVAR) w.r.t. the j-th parameter, j=1...NTLM KPP_REAL :: Y_tlm(NVAR,NTLM) KPP_REAL :: TIN ! TIN - Start Time KPP_REAL :: TOUT ! TOUT - End Time KPP_REAL, INTENT(IN), OPTIONAL :: RTOL_tlm(NVAR,NTLM),ATOL_tlm(NVAR,NTLM) INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) INTEGER, INTENT(OUT), OPTIONAL :: IERR_U INTEGER :: IERR KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2 INTEGER :: ICNTRL(20), ISTATUS(20) INTEGER, SAVE :: Ntotal = 0 ICNTRL(1:20) = 0 RCNTRL(1:20) = 0.0_dp !~~~> fine-tune the integrator: ICNTRL(2) = 0 ! 0=vector tolerances, 1=scalar tolerances ICNTRL(5) = 8 ! Max no. of Newton iterations ICNTRL(6) = 0 ! Starting values for Newton are interpolated (0) or zero (1) ICNTRL(7) = 0 ! How to solve TLM: 0=modified Newton, 1=direct ICNTRL(9) = 0 ! TLM Newton Iterations influence ICNTRL(12) = 0 ! TLM Truncation Error influence !~~~> if optional parameters are given, and if they are >0, ! then use them to overwrite default settings IF (PRESENT(ICNTRL_U)) THEN WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:) END IF IF (PRESENT(RCNTRL_U)) THEN WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:) END IF T1 = TIN; T2 = TOUT CALL SdirkTLM( NVAR, NTLM, T1, T2, Y, Y_tlm, RTOL, ATOL, & RTOL_tlm, ATOL_tlm, RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR ) !~~~> Debug option: print number of steps ! Ntotal = Ntotal + ISTATUS(Nstp) ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')',' O3=', VAR(ind_O3) ! if optional parameters are given for output ! use them to store information in them IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:) IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:) IF (PRESENT(IERR_U)) IERR_U = IERR IF (IERR < 0) THEN PRINT *,'SDIRK-TLM: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')' ENDIF END SUBROUTINE INTEGRATE_TLM !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SdirkTLM(N, NTLM, Tinitial, Tfinal, Y, Y_tlm, RelTol, AbsTol, & RelTol_tlm, AbsTol_tlm, & RCNTRL, ICNTRL, RSTATUS, ISTATUS, Ierr) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit ! Runge-Kutta (SDIRK) method. ! ! This implementation is based on the book and the code Sdirk4: ! ! E. Hairer and G. Wanner ! "Solving ODEs II. Stiff and differential-algebraic problems". ! Springer series in computational mathematics, Springer-Verlag, 1996. ! This code is based on the SDIRK4 routine in the above book. ! ! Methods: ! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 ! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant ! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 ! ! (C) Adrian Sandu, July 2005 ! Virginia Polytechnic Institute and State University ! Contact: sandu@cs.vt.edu ! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! This implementation is part of KPP - the Kinetic PreProcessor !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> INPUT ARGUMENTS: ! !- Y(NVAR) = vector of initial conditions (at T=Tinitial) !- [Tinitial,Tfinal] = time range of integration ! (if Tinitial>Tfinal the integration is performed backwards in time) !- RelTol, AbsTol = user precribed accuracy !- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function, ! returns Ydot = Y' = F(T,Y) !- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function, ! returns Jcb = dF/dY !- ICNTRL(1:20) = integer inputs parameters !- RCNTRL(1:20) = real inputs parameters !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> OUTPUT ARGUMENTS: ! !- Y(NVAR) -> vector of final states (at T->Tfinal) !- ISTATUS(1:20) -> integer output parameters !- RSTATUS(1:20) -> real output parameters !- Ierr -> job status upon return ! success (positive value) or ! failure (negative value) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> INPUT PARAMETERS: ! ! Note: For input parameters equal to zero the default values of the ! corresponding variables are used. ! ! Note: For input parameters equal to zero the default values of the ! corresponding variables are used. !~~~> ! ICNTRL(1) = not used ! ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors ! = 1: AbsTol, RelTol are scalars ! ! ICNTRL(3) = Method ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0 the default value of 100000 is used ! ! ICNTRL(5) -> maximum number of Newton iterations ! For ICNTRL(5)=0 the default value of 8 is used ! ! ICNTRL(6) -> starting values of Newton iterations: ! ICNTRL(6)=0 : starting values are interpolated (the default) ! ICNTRL(6)=1 : starting values are zero ! ! ICNTRL(7) -> method to solve TLM equations: ! ICNTRL(7)=0 : modified Newton re-using LU (the default) ! ICNTRL(7)=1 : direct solution (additional one LU factorization per stage) ! ! ICNTRL(9) -> switch for TLM Newton iteration error estimation strategy ! ICNTRL(9) = 0: base number of iterations as forward solution ! ICNTRL(9) = 1: use RTOL_tlm and ATOL_tlm to calculate ! error estimation for TLM at Newton stages ! ! ICNTRL(12) -> switch for TLM truncation error control ! ICNTRL(12) = 0: TLM error is not used ! ICNTRL(12) = 1: TLM error is computed and used ! ! !~~~> Real parameters ! ! RCNTRL(1) -> Hmin, lower bound for the integration step size ! It is strongly recommended to keep Hmin = ZERO ! RCNTRL(2) -> Hmax, upper bound for the integration step size ! RCNTRL(3) -> Hstart, starting value for the integration step size ! ! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2) ! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6) ! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections ! (default=0.1) ! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller ! than the predicted value (default=0.9) ! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller ! than ThetaMin the Jacobian is not recomputed; ! (default=0.001) ! RCNTRL(9) -> NewtonTol, stopping criterion for Newton's method ! (default=0.03) ! RCNTRL(10) -> Qmin ! RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the ! step size is kept constant and the LU factorization ! reused (default Qmin=1, Qmax=1.2) ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> OUTPUT PARAMETERS: ! ! Note: each call to Rosenbrock adds the current no. of fcn calls ! to previous value of ISTATUS(1), and similar for the other params. ! Set ISTATUS(1:10) = 0 before call to avoid this accumulation. ! ! ISTATUS(1) = No. of function calls ! ISTATUS(2) = No. of jacobian calls ! ISTATUS(3) = No. of steps ! ISTATUS(4) = No. of accepted steps ! ISTATUS(5) = No. of rejected steps (except at the beginning) ! ISTATUS(6) = No. of LU decompositions ! ISTATUS(7) = No. of forward/backward substitutions ! ISTATUS(8) = No. of singular matrix decompositions ! ! RSTATUS(1) -> Texit, the time corresponding to the ! computed Y upon return ! RSTATUS(2) -> Hexit,last accepted step before return ! RSTATUS(3) -> Hnew, last predicted step before return ! For multiple restarts, use Hnew as Hstart in the following run ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE ! Arguments INTEGER, INTENT(IN) :: N, NTLM, ICNTRL(20) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, & RelTol(N), AbsTol(N), RCNTRL(20), & RelTol_tlm(N,NTLM), AbsTol_tlm(N,NTLM) KPP_REAL, INTENT(INOUT) :: Y(NVAR), Y_tlm(N,NTLM) INTEGER, INTENT(OUT) :: Ierr INTEGER, INTENT(INOUT) :: ISTATUS(20) KPP_REAL, INTENT(OUT) :: RSTATUS(20) !~~~> SDIRK method coefficients, up to 5 stages INTEGER, PARAMETER :: Smax = 5 INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5 KPP_REAL :: rkGamma, rkA(Smax,Smax), rkB(Smax), rkC(Smax), & rkD(Smax), rkE(Smax), rkBhat(Smax), rkELO, & rkAlpha(Smax,Smax), rkTheta(Smax,Smax) INTEGER :: sdMethod, rkS ! The number of stages ! Local variables KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, & FacMin, Facmax, FacSafe, FacRej, & ThetaMin, NewtonTol, Qmin, Qmax INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i LOGICAL :: StartNewton, DirectTLM, TLMNewtonEst, TLMtruncErr KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 Ierr = 0 ISTATUS(1:20) = 0 RSTATUS(1:20) = ZERO !~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1) ! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR) IF (ICNTRL(2) == 0) THEN ITOL = 1 ELSE ITOL = 0 END IF !~~~> ICNTRL(3) - method selection SELECT CASE (ICNTRL(3)) CASE (0,1) CALL Sdirk2a CASE (2) CALL Sdirk2b CASE (3) CALL Sdirk3a CASE (4) CALL Sdirk4a CASE (5) CALL Sdirk4b CASE DEFAULT CALL Sdirk2a END SELECT !~~~> The maximum number of time steps admitted IF (ICNTRL(4) == 0) THEN Max_no_steps = 200000 ELSEIF (ICNTRL(4) > 0) THEN Max_no_steps=ICNTRL(4) ELSE PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4) CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr) END IF !~~~> The maximum number of Newton iterations admitted IF(ICNTRL(5) == 0)THEN NewtonMaxit=8 ELSE NewtonMaxit=ICNTRL(5) IF(NewtonMaxit <= 0)THEN PRINT * ,'User-selected ICNTRL(5)=',ICNTRL(5) CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,Ierr) END IF END IF !~~~> StartNewton: Use extrapolation for starting values of Newton iterations IF (ICNTRL(6) == 0) THEN StartNewton = .TRUE. ELSE StartNewton = .FALSE. END IF !~~~> Solve TLM equations directly or by Newton iterations DirectTLM = (ICNTRL(7) == 1) !~~~> Newton iteration error control selection IF (ICNTRL(9) == 0) THEN TLMNewtonEst = .FALSE. ELSE TLMNewtonEst = .TRUE. END IF !~~~> TLM truncation error control selection IF (ICNTRL(12) == 0) THEN TLMtruncErr = .FALSE. ELSE TLMtruncErr = .TRUE. END IF !~~~> Unit roundoff (1+Roundoff>1) Roundoff = WLAMCH('E') !~~~> Lower bound on the step size: (positive value) IF (RCNTRL(1) == ZERO) THEN Hmin = ZERO ELSEIF (RCNTRL(1) > ZERO) THEN Hmin = RCNTRL(1) ELSE PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF !~~~> Upper bound on the step size: (positive value) IF (RCNTRL(2) == ZERO) THEN Hmax = ABS(Tfinal-Tinitial) ELSEIF (RCNTRL(2) > ZERO) THEN Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial)) ELSE PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN Hstart = MAX(Hmin,Roundoff) ELSEIF (RCNTRL(3) > ZERO) THEN Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial)) ELSE PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr) END IF !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax IF (RCNTRL(4) == ZERO) THEN FacMin = 0.2_dp ELSEIF (RCNTRL(4) > ZERO) THEN FacMin = RCNTRL(4) ELSE PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4) CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) END IF IF (RCNTRL(5) == ZERO) THEN FacMax = 10.0_dp ELSEIF (RCNTRL(5) > ZERO) THEN FacMax = RCNTRL(5) ELSE PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5) CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) END IF !~~~> FacRej: Factor to decrease step after 2 succesive rejections IF (RCNTRL(6) == ZERO) THEN FacRej = 0.1_dp ELSEIF (RCNTRL(6) > ZERO) THEN FacRej = RCNTRL(6) ELSE PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6) CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) END IF !~~~> FacSafe: Safety Factor in the computation of new step size IF (RCNTRL(7) == ZERO) THEN FacSafe = 0.9_dp ELSEIF (RCNTRL(7) > ZERO) THEN FacSafe = RCNTRL(7) ELSE PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7) CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr) END IF !~~~> ThetaMin: decides whether the Jacobian should be recomputed IF(RCNTRL(8) == 0.D0)THEN ThetaMin = 1.0d-3 ELSE ThetaMin = RCNTRL(8) END IF !~~~> Stopping criterion for Newton's method IF(RCNTRL(9) == ZERO)THEN NewtonTol = 3.0d-2 ELSE NewtonTol = RCNTRL(9) END IF !~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST. IF(RCNTRL(10) == ZERO)THEN Qmin=ONE ELSE Qmin=RCNTRL(10) END IF IF(RCNTRL(11) == ZERO)THEN Qmax=1.2D0 ELSE Qmax=RCNTRL(11) END IF !~~~> Check if tolerances are reasonable IF (ITOL == 0) THEN IF (AbsTol(1) <= ZERO .OR. RelTol(1) <= 10.D0*Roundoff) THEN PRINT * , ' Scalar AbsTol = ',AbsTol(1) PRINT * , ' Scalar RelTol = ',RelTol(1) CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr) END IF ELSE DO i=1,N IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN PRINT * , ' AbsTol(',i,') = ',AbsTol(i) PRINT * , ' RelTol(',i,') = ',RelTol(i) CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr) END IF END DO END IF IF (Ierr < 0) RETURN CALL SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CONTAINS ! PROCEDURES INTERNAL TO SDIRK !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE KPP_ROOT_Parameters IMPLICIT NONE !~~~> Arguments: INTEGER, INTENT(IN) :: N, NTLM KPP_REAL, INTENT(INOUT) :: Y(N), Y_tlm(N,NTLM) KPP_REAL, INTENT(IN) :: Tinitial, Tfinal INTEGER, INTENT(OUT) :: Ierr !~~~> Local variables: KPP_REAL :: Z(NVAR,rkS), G(NVAR), TMP(NVAR), & NewtonRate, SCAL(NVAR), DZ(NVAR), & T, H, Theta, Hratio, NewtonPredictedErr, & Qnewton, Err, Fac, Hnew, Tdirection, & NewtonIncrement, NewtonIncrementOld, & SCAL_tlm(NVAR), Yerr(N), Yerr_tlm(N,NTLM), ThetaTLM KPP_REAL :: Z_tlm(NVAR,rkS,NTLM) INTEGER :: itlm, j, IER, istage, NewtonIter, saveNiter, NewtonIterTLM INTEGER :: IP(NVAR), IP_tlm(NVAR) LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone #ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: FJAC, E, Jac, E_tlm #else KPP_REAL, DIMENSION(LU_NONZERO) :: FJAC, E, Jac, E_tlm #endif KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Initializations !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ T = Tinitial Tdirection = SIGN(ONE,Tfinal-Tinitial) H = MAX(ABS(Hmin),ABS(Hstart)) IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6 H=MIN(ABS(H),Hmax) H=SIGN(H,Tdirection) SkipLU = .FALSE. SkipJac = .FALSE. Reject = .FALSE. FirstStep=.TRUE. CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Time loop begins !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO ) !~~~> Compute E = 1/(h*gamma)-Jac and its LU decomposition IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipJac, SkipLU, E, IP, Reject, IER ) IF (IER /= 0) THEN CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN END IF END IF IF (ISTATUS(Nstp) > Max_no_steps) THEN CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN END IF IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN END IF stages:DO istage = 1, rkS !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Simplified Newton iterations !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Starting values for Newton iterations CALL Set2zero(N,Z(1,istage)) !~~~> Prepare the loop-independent part of the right-hand side CALL Set2zero(N,G) IF (istage > 1) THEN DO j = 1, istage-1 ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:)) CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1) ! Zi(:) = sum_j Alpha(i,j)*Zj(:) IF (StartNewton) THEN CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1) END IF END DO END IF !~~~> Initializations for Newton iteration NewtonDone = .FALSE. Fac = 0.5d0 ! Step reduction factor if too many iterations NewtonLoop:DO NewtonIter = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi CALL FUN_CHEM(T+rkC(istage)*H,TMP,DZ) ! DZ <- Fun(Y+Zi) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 ! DZ(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*DZ(1:N) CALL WSCAL(N, H*rkGamma, DZ, 1) CALL WAXPY (N, -ONE, Z(1,istage), 1, DZ, 1) CALL WAXPY (N, ONE, G,1, DZ,1) !~~~> Solve the linear system CALL SDIRK_Solve ( H, N, E, IP, IER, DZ ) !~~~> Check convergence of Newton iterations CALL SDIRK_ErrorNorm(N, DZ, SCAL, NewtonIncrement) IF ( NewtonIter == 1 ) THEN Theta = ABS(ThetaMin) NewtonRate = 2.0d0 ELSE Theta = NewtonIncrement/NewtonIncrementOld IF (Theta < 0.99d0) THEN NewtonRate = Theta/(ONE-Theta) ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta) IF (NewtonPredictedErr >= NewtonTol) THEN ! Non-convergence of Newton: predicted error too large Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol) Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter)) EXIT NewtonLoop END IF ELSE ! Non-convergence of Newton: Theta too large EXIT NewtonLoop END IF END IF NewtonIncrementOld = NewtonIncrement ! Update solution: Z(:) <-- Z(:)+DZ(:) CALL WAXPY(N,ONE,DZ,1,Z(1,istage),1) ! Check error in Newton iterations NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) THEN ! Tune error in TLM variables by defining the minimal number of Newton iterations. saveNiter = NewtonIter+1 EXIT NewtonLoop END IF END DO NewtonLoop IF (.NOT.NewtonDone) THEN !CALL RK_ErrorMsg(-12,T,H,Ierr); H = Fac*H; Reject=.TRUE. SkipJac = .TRUE.; SkipLU = .FALSE. CYCLE Tloop END IF !~~~> End of simplified Newton iterations for forward variables !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Solve for TLM variables !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Direct solution for TLM variables DirTLM:IF (DirectTLM) THEN TMP(1:N) = Y(1:N) + Z(1:N,istage) SkipJac = .FALSE. CALL SDIRK_PrepareMatrix ( H, T+rkC(istage)*H, TMP, Jac, & SkipJac, SkipLU, E_tlm, IP_tlm, Reject, IER ) IF (IER /= 0) CYCLE TLoop TlmL: DO itlm = 1, NTLM G(1:N) = Y_tlm(1:N,itlm) IF (istage > 1) THEN ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:) ! = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm) DO j = 1, istage-1 CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1) END DO END IF CALL SDIRK_Solve ( H, N, E_tlm, IP_tlm, IER, G ) Z_tlm(1:N,istage,itlm) = G(1:N) - Y_tlm(1:N,itlm) END DO TlmL ELSE DirTLM !~~~> Jacobian of the current stage solution TMP(1:N) = Y(1:N) + Z(1:N,istage) CALL JAC_CHEM(T+rkC(istage)*H,TMP,Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 !~~~> Simplified Newton iterations for TLM variables TlmLoop:DO itlm = 1,NTLM NewtonRate = MAX(NewtonRate,Roundoff)**0.8d0 !~~~> Starting values for Newton iterations CALL Set2zero(N,Z_tlm(1,istage,itlm)) !~~~> Prepare the loop-independent part of the right-hand side #ifdef FULL_ALGEBRA DZ = MATMUL(Jac,Y_tlm(1,itlm)) ! DZ <- Jac(Y+Z)*Y_tlm #else CALL Jac_SP_Vec ( Jac, Y_tlm(1,itlm), DZ ) #endif G(1:N) = (H*rkGamma)*DZ(1:N) IF (istage > 1) THEN ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:) ! = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm) DO j = 1, istage-1 CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1) END DO END IF !~~~> Initializations for Newton iteration IF (TLMNewtonEst) THEN NewtonDone = .FALSE. Fac = 0.5d0 ! Step reduction factor if too many iterations CALL SDIRK_ErrorScale(N,ITOL,AbsTol_tlm(1,itlm),RelTol_tlm(1,itlm), & Y_tlm(1,itlm),SCAL_tlm) END IF NewtonLoopTLM:DO NewtonIterTLM = 1, NewtonMaxit !~~~> Prepare the loop-dependent part of the right-hand side #ifdef FULL_ALGEBRA DZ = MATMUL(Jac,Z_tlm(1,istage,itlm)) ! DZ <- Jac(Y+Z)*Z_tlm #else CALL Jac_SP_Vec ( Jac, Z_tlm(1,istage,itlm), DZ ) #endif DZ(1:N) = (H*rkGamma)*DZ(1:N)+G(1:N)-Z_tlm(1:N,istage,itlm) CALL SDIRK_Solve ( H, N, E, IP, IER, DZ ) IF (TLMNewtonEst) THEN !~~~> Check convergence of Newton iterations CALL SDIRK_ErrorNorm(N, DZ, SCAL_tlm, NewtonIncrement) IF ( NewtonIterTLM <= 1 ) THEN ThetaTLM = ABS(ThetaMin) NewtonRate = 2.0d0 ELSE ThetaTLM = NewtonIncrement/NewtonIncrementOld IF (ThetaTLM < 0.99d0) THEN NewtonRate = ThetaTLM/(ONE-ThetaTLM) ! Predict error at the end of Newton process NewtonPredictedErr = NewtonIncrement & *ThetaTLM**(NewtonMaxit-NewtonIterTLM)/(ONE-ThetaTLM) IF (NewtonPredictedErr >= NewtonTol) THEN ! Non-convergence of Newton: predicted error too large Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol) Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIterTLM)) EXIT NewtonLoopTLM END IF ELSE ! Non-convergence of Newton: Theta too large EXIT NewtonLoopTLM END IF END IF NewtonIncrementOld = NewtonIncrement END IF !(TLMNewtonEst) ! Update solution: Z_tlm(:) <-- Z_tlm(:)+DZ(:) CALL WAXPY(N,ONE,DZ,1,Z_tlm(1,istage,itlm),1) ! Check error in Newton iterations IF (TLMNewtonEst) THEN NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol) IF (NewtonDone) EXIT NewtonLoopTLM ELSE ! Minimum number of iterations same as FWD iterations IF (NewtonIterTLM>=saveNiter) EXIT NewtonLoopTLM END IF END DO NewtonLoopTLM IF ((TLMNewtonEst) .AND. (.NOT.NewtonDone)) THEN !CALL RK_ErrorMsg(-12,T,H,Ierr); H = Fac*H; Reject=.TRUE. SkipJac = .TRUE.; SkipLU = .FALSE. CYCLE Tloop END IF END DO TlmLoop !~~~> End of simplified Newton iterations for TLM END IF DirTLM END DO stages !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Error estimation !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ISTATUS(Nstp) = ISTATUS(Nstp) + 1 CALL Set2zero(N,Yerr) DO i = 1,rkS IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,Yerr,1) END DO CALL SDIRK_Solve ( H, N, E, IP, IER, Yerr ) CALL SDIRK_ErrorNorm(N, Yerr, SCAL, Err) IF (TLMtruncErr) THEN CALL Set2zero(NVAR*NTLM,Yerr_tlm) DO itlm=1,NTLM DO j=1,rkS IF (rkE(j) /= ZERO) CALL WAXPY(N,rkE(j),Z_tlm(1,j,itlm),1,Yerr_tlm(1,itlm),1) END DO CALL SDIRK_Solve (H, N, E, IP, IER, Yerr_tlm(1,itlm)) END DO CALL SDIRK_ErrorNorm_tlm(N,NTLM, Yerr_tlm, Err) END IF !~~~> Computation of new step size Hnew Fac = FacSafe*(Err)**(-ONE/rkELO) Fac = MAX(FacMin,MIN(FacMax,Fac)) Hnew = H*Fac !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Accept/Reject step !~~~~~~~~~~~~~~~~~~~~~~~~~~~~ accept: IF ( Err < ONE ) THEN !~~~> Step is accepted FirstStep=.FALSE. ISTATUS(Nacc) = ISTATUS(Nacc) + 1 !~~~> Update time and solution T = T + H ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:) DO i = 1,rkS IF (rkD(i)/=ZERO) THEN CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1) DO itlm = 1, NTLM CALL WAXPY(N,rkD(i),Z_tlm(1,i,itlm),1,Y_tlm(1,itlm),1) END DO END IF END DO !~~~> Update scaling coefficients CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) !~~~> Next time step Hnew = Tdirection*MIN(ABS(Hnew),Hmax) ! Last T and H RSTATUS(Ntexit) = T RSTATUS(Nhexit) = H RSTATUS(Nhnew) = Hnew ! No step increase after a rejection IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H)) Reject = .FALSE. IF ((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO) THEN H = Tfinal-T ELSE Hratio=Hnew/H ! If step not changed too much keep Jacobian and reuse LU SkipLU = ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) & .AND. (Hratio <= Qmax) ) ! For TLM: do not skip LU (decrease TLM error) SkipLU = .FALSE. IF (.NOT.SkipLU) H = Hnew END IF ! If convergence is fast enough, do not update Jacobian ! SkipJac = (Theta <= ThetaMin) SkipJac = .FALSE. ELSE accept !~~~> Step is rejected IF (FirstStep .OR. Reject) THEN H = FacRej*H ELSE H = Hnew END IF Reject = .TRUE. SkipJac = .TRUE. SkipLU = .FALSE. IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1 END IF accept END DO Tloop ! Successful return Ierr = 1 END SUBROUTINE SDIRK_IntegratorTLM !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE INTEGER :: N, i, ITOL KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), & Y(NVAR), SCAL(NVAR) IF (ITOL == 0) THEN DO i=1,N SCAL(i) = ONE / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) ) END DO ELSE DO i=1,N SCAL(i) = ONE / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) ) END DO END IF END SUBROUTINE SDIRK_ErrorScale !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! INTEGER :: i, N KPP_REAL :: Y(N), SCAL(N), Err Err = ZERO DO i=1,N Err = Err+(Y(i)*SCAL(i))**2 END DO Err = MAX( SQRT(Err/DBLE(N)), 1.0d-10 ) ! END SUBROUTINE SDIRK_ErrorNorm !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorNorm_tlm(N,NTLM, Y_tlm, FWD_Err) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! INTEGER :: itlm, NTLM, N KPP_REAL :: Y_tlm(N,NTLM), SCAL_tlm(N), FWD_Err, Err DO itlm=1,NTLM CALL SDIRK_ErrorScale(N,ITOL,AbsTol_tlm(1,itlm),RelTol_tlm(1,itlm), & Y_tlm(1,itlm),SCAL_tlm) CALL SDIRK_ErrorNorm(N, Y_tlm(1,itlm), SCAL_tlm, Err) FWD_Err = MAX(FWD_Err, Err) END DO ! END SUBROUTINE SDIRK_ErrorNorm_tlm !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Handles all error messages !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ KPP_REAL, INTENT(IN) :: T, H INTEGER, INTENT(IN) :: Code INTEGER, INTENT(OUT) :: Ierr Ierr = Code PRINT * , & 'Forced exit from SDIRK due to the following error:' SELECT CASE (Code) CASE (-1) PRINT * , '--> Improper value for maximal no of steps' CASE (-2) PRINT * , '--> Improper value for maximal no of Newton iterations' CASE (-3) PRINT * , '--> Hmin/Hmax/Hstart must be positive' CASE (-4) PRINT * , '--> FacMin/FacMax/FacRej must be positive' CASE (-5) PRINT * , '--> Improper tolerance values' CASE (-6) PRINT * , '--> No of steps exceeds maximum bound' CASE (-7) PRINT * , '--> Step size too small: T + 10*H = T', & ' or H < Roundoff' CASE (-8) PRINT * , '--> Matrix is repeatedly singular' CASE DEFAULT PRINT *, 'Unknown Error code: ', Code END SELECT PRINT *, "T=", T, "and H=", H END SUBROUTINE SDIRK_ErrorMsg !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, & SkipJac, SkipLU, E, IP, Reject, ISING ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Compute the matrix E = I - 1/(H*Gamma)*Jac, and its decomposition !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE KPP_REAL, INTENT(INOUT) :: H KPP_REAL, INTENT(IN) :: T, Y(NVAR) LOGICAL, INTENT(INOUT) :: SkipJac,SkipLU,Reject INTEGER, INTENT(OUT) :: ISING, IP(NVAR) #ifdef FULL_ALGEBRA KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR) KPP_REAL, INTENT(OUT) :: E(NVAR,NVAR) #else KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO) KPP_REAL, INTENT(OUT) :: E(LU_NONZERO) #endif KPP_REAL :: HGammaInv INTEGER :: i, j, ConsecutiveSng ConsecutiveSng = 0 ISING = 1 Hloop: DO WHILE (ISING /= 0) HGammaInv = ONE/(H*rkGamma) !~~~> Compute the Jacobian ! IF (SkipJac) THEN ! SkipJac = .FALSE. ! ELSE IF (.NOT. SkipJac) THEN CALL JAC_CHEM( T, Y, FJAC ) ISTATUS(Njac) = ISTATUS(Njac) + 1 END IF #ifdef FULL_ALGEBRA DO j=1,NVAR DO i=1,NVAR E(i,j) = -FJAC(i,j) END DO E(j,j) = E(j,j)+HGammaInv END DO CALL DGETRF( NVAR, NVAR, E, NVAR, IP, ISING ) #else DO i = 1,LU_NONZERO E(i) = -FJAC(i) END DO DO i = 1,NVAR j = LU_DIAG(i); E(j) = E(j) + HGammaInv END DO CALL KppDecomp ( E, ISING) IP(1) = 1 #endif ISTATUS(Ndec) = ISTATUS(Ndec) + 1 IF (ISING /= 0) THEN WRITE (6,*) ' MATRIX IS SINGULAR, ISING=',ISING,' T=',T,' H=',H ISTATUS(Nsng) = ISTATUS(Nsng) + 1; ConsecutiveSng = ConsecutiveSng + 1 IF (ConsecutiveSng >= 6) RETURN ! Failure H = 0.5d0*H SkipJac = .TRUE. SkipLU = .FALSE. Reject = .TRUE. END IF END DO Hloop END SUBROUTINE SDIRK_PrepareMatrix !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Solves the system (H*Gamma-Jac)*x = R ! using the LU decomposition of E = I - 1/(H*Gamma)*Jac !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE INTEGER, INTENT(IN) :: N, IP(N), ISING KPP_REAL, INTENT(IN) :: H #ifdef FULL_ALGEBRA KPP_REAL, INTENT(IN) :: E(NVAR,NVAR) #else KPP_REAL, INTENT(IN) :: E(LU_NONZERO) #endif KPP_REAL, INTENT(INOUT) :: RHS(N) KPP_REAL :: HGammaInv HGammaInv = ONE/(H*rkGamma) CALL WSCAL(N,HGammaInv,RHS,1) #ifdef FULL_ALGEBRA CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING ) #else CALL KppSolve(E, RHS) #endif ISTATUS(Nsol) = ISTATUS(Nsol) + 1 END SUBROUTINE SDIRK_Solve !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Sdirk4a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S4A ! Number of stages rkS = 5 ! Method coefficients rkGamma = .2666666666666666666666666666666667d0 rkA(1,1) = .2666666666666666666666666666666667d0 rkA(2,1) = .5000000000000000000000000000000000d0 rkA(2,2) = .2666666666666666666666666666666667d0 rkA(3,1) = .3541539528432732316227461858529820d0 rkA(3,2) = -.5415395284327323162274618585298197d-1 rkA(3,3) = .2666666666666666666666666666666667d0 rkA(4,1) = .8515494131138652076337791881433756d-1 rkA(4,2) = -.6484332287891555171683963466229754d-1 rkA(4,3) = .7915325296404206392428857585141242d-1 rkA(4,4) = .2666666666666666666666666666666667d0 rkA(5,1) = 2.100115700566932777970612055999074d0 rkA(5,2) = -.7677800284445976813343102185062276d0 rkA(5,3) = 2.399816361080026398094746205273880d0 rkA(5,4) = -2.998818699869028161397714709433394d0 rkA(5,5) = .2666666666666666666666666666666667d0 rkB(1) = 2.100115700566932777970612055999074d0 rkB(2) = -.7677800284445976813343102185062276d0 rkB(3) = 2.399816361080026398094746205273880d0 rkB(4) = -2.998818699869028161397714709433394d0 rkB(5) = .2666666666666666666666666666666667d0 rkBhat(1)= 2.885264204387193942183851612883390d0 rkBhat(2)= -.1458793482962771337341223443218041d0 rkBhat(3)= 2.390008682465139866479830743628554d0 rkBhat(4)= -4.129393538556056674929560012190140d0 rkBhat(5)= 0.d0 rkC(1) = .2666666666666666666666666666666667d0 rkC(2) = .7666666666666666666666666666666667d0 rkC(3) = .5666666666666666666666666666666667d0 rkC(4) = .3661315380631796996374935266701191d0 rkC(5) = 1.d0 ! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} rkD(1) = 0.d0 rkD(2) = 0.d0 rkD(3) = 0.d0 rkD(4) = 0.d0 rkD(5) = 1.d0 ! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} rkE(1) = -.6804000050475287124787034884002302d0 rkE(2) = 1.558961944525217193393931795738823d0 rkE(3) = -13.55893003128907927748632408763868d0 rkE(4) = 15.48522576958521253098585004571302d0 rkE(5) = 1.d0 ! Local order of Err estimate rkElo = 4 ! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} rkTheta(2,1) = 1.875000000000000000000000000000000d0 rkTheta(3,1) = 1.708847304091539528432732316227462d0 rkTheta(3,2) = -.2030773231622746185852981969486824d0 rkTheta(4,1) = .2680325578937783958847157206823118d0 rkTheta(4,2) = -.1828840955527181631794050728644549d0 rkTheta(4,3) = .2968246986151577397160821594427966d0 rkTheta(5,1) = .9096171815241460655379433581446771d0 rkTheta(5,2) = -3.108254967778352416114774430509465d0 rkTheta(5,3) = 12.33727431701306195581826123274001d0 rkTheta(5,4) = -11.24557012450885560524143016037523d0 ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 2.875000000000000000000000000000000d0 rkAlpha(3,1) = .8500000000000000000000000000000000d0 rkAlpha(3,2) = .4434782608695652173913043478260870d0 rkAlpha(4,1) = .7352046091658870564637910527807370d0 rkAlpha(4,2) = -.9525565003057343527941920657462074d-1 rkAlpha(4,3) = .4290111305453813852259481840631738d0 rkAlpha(5,1) = -16.10898993405067684831655675112808d0 rkAlpha(5,2) = 6.559571569643355712998131800797873d0 rkAlpha(5,3) = -15.90772144271326504260996815012482d0 rkAlpha(5,4) = 25.34908987169226073668861694892683d0 !~~~> Coefficients for continuous solution ! rkD(1,1)= 24.74416644927758d0 ! rkD(1,2)= -4.325375951824688d0 ! rkD(1,3)= 41.39683763286316d0 ! rkD(1,4)= -61.04144619901784d0 ! rkD(1,5)= -3.391332232917013d0 ! rkD(2,1)= -51.98245719616925d0 ! rkD(2,2)= 10.52501981094525d0 ! rkD(2,3)= -154.2067922191855d0 ! rkD(2,4)= 214.3082125319825d0 ! rkD(2,5)= 14.71166018088679d0 ! rkD(3,1)= 33.14347947522142d0 ! rkD(3,2)= -19.72986789558523d0 ! rkD(3,3)= 230.4878502285804d0 ! rkD(3,4)= -287.6629744338197d0 ! rkD(3,5)= -18.99932366302254d0 ! rkD(4,1)= -5.905188728329743d0 ! rkD(4,2)= 13.53022403646467d0 ! rkD(4,3)= -117.6778956422581d0 ! rkD(4,4)= 134.3962081008550d0 ! rkD(4,5)= 8.678995715052762d0 ! ! DO i=1,4 ! CONTi <-- Sum_j rkD(i,j)*Zj ! CALL Set2zero(N,CONT(1,i)) ! DO j = 1,rkS ! CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1) ! END DO ! END DO rkELO = 4.0d0 END SUBROUTINE Sdirk4a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Sdirk4b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S4B ! Number of stages rkS = 5 ! Method coefficients rkGamma = .25d0 rkA(1,1) = 0.25d0 rkA(2,1) = 0.5d00 rkA(2,2) = 0.25d0 rkA(3,1) = 0.34d0 rkA(3,2) =-0.40d-1 rkA(3,3) = 0.25d0 rkA(4,1) = 0.2727941176470588235294117647058824d0 rkA(4,2) =-0.5036764705882352941176470588235294d-1 rkA(4,3) = 0.2757352941176470588235294117647059d-1 rkA(4,4) = 0.25d0 rkA(5,1) = 1.041666666666666666666666666666667d0 rkA(5,2) =-1.020833333333333333333333333333333d0 rkA(5,3) = 7.812500000000000000000000000000000d0 rkA(5,4) =-7.083333333333333333333333333333333d0 rkA(5,5) = 0.25d0 rkB(1) = 1.041666666666666666666666666666667d0 rkB(2) = -1.020833333333333333333333333333333d0 rkB(3) = 7.812500000000000000000000000000000d0 rkB(4) = -7.083333333333333333333333333333333d0 rkB(5) = 0.250000000000000000000000000000000d0 rkBhat(1)= 1.069791666666666666666666666666667d0 rkBhat(2)= -0.894270833333333333333333333333333d0 rkBhat(3)= 7.695312500000000000000000000000000d0 rkBhat(4)= -7.083333333333333333333333333333333d0 rkBhat(5)= 0.212500000000000000000000000000000d0 rkC(1) = 0.25d0 rkC(2) = 0.75d0 rkC(3) = 0.55d0 rkC(4) = 0.50d0 rkC(5) = 1.00d0 ! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} rkD(1) = 0.0d0 rkD(2) = 0.0d0 rkD(3) = 0.0d0 rkD(4) = 0.0d0 rkD(5) = 1.0d0 ! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} rkE(1) = 0.5750d0 rkE(2) = 0.2125d0 rkE(3) = -4.6875d0 rkE(4) = 4.2500d0 rkE(5) = 0.1500d0 ! Local order of Err estimate rkElo = 4 ! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} rkTheta(2,1) = 2.d0 rkTheta(3,1) = 1.680000000000000000000000000000000d0 rkTheta(3,2) = -.1600000000000000000000000000000000d0 rkTheta(4,1) = 1.308823529411764705882352941176471d0 rkTheta(4,2) = -.1838235294117647058823529411764706d0 rkTheta(4,3) = 0.1102941176470588235294117647058824d0 rkTheta(5,1) = -3.083333333333333333333333333333333d0 rkTheta(5,2) = -4.291666666666666666666666666666667d0 rkTheta(5,3) = 34.37500000000000000000000000000000d0 rkTheta(5,4) = -28.33333333333333333333333333333333d0 ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 3. rkAlpha(3,1) = .8800000000000000000000000000000000d0 rkAlpha(3,2) = .4400000000000000000000000000000000d0 rkAlpha(4,1) = .1666666666666666666666666666666667d0 rkAlpha(4,2) = -.8333333333333333333333333333333333d-1 rkAlpha(4,3) = .9469696969696969696969696969696970d0 rkAlpha(5,1) = -6.d0 rkAlpha(5,2) = 9.d0 rkAlpha(5,3) = -56.81818181818181818181818181818182d0 rkAlpha(5,4) = 54.d0 END SUBROUTINE Sdirk4b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Sdirk2a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S2A ! Number of stages rkS = 2 ! Method coefficients rkGamma = .2928932188134524755991556378951510d0 rkA(1,1) = .2928932188134524755991556378951510d0 rkA(2,1) = .7071067811865475244008443621048490d0 rkA(2,2) = .2928932188134524755991556378951510d0 rkB(1) = .7071067811865475244008443621048490d0 rkB(2) = .2928932188134524755991556378951510d0 rkBhat(1)= .6666666666666666666666666666666667d0 rkBhat(2)= .3333333333333333333333333333333333d0 rkC(1) = 0.292893218813452475599155637895151d0 rkC(2) = 1.0d0 ! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} rkD(1) = 0.0d0 rkD(2) = 1.0d0 ! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} rkE(1) = 0.4714045207910316829338962414032326d0 rkE(2) = -0.1380711874576983496005629080698993d0 ! Local order of Err estimate rkElo = 2 ! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} rkTheta(2,1) = 2.414213562373095048801688724209698d0 ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 3.414213562373095048801688724209698d0 END SUBROUTINE Sdirk2a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Sdirk2b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S2B ! Number of stages rkS = 2 ! Method coefficients rkGamma = 1.707106781186547524400844362104849d0 rkA(1,1) = 1.707106781186547524400844362104849d0 rkA(2,1) = -.707106781186547524400844362104849d0 rkA(2,2) = 1.707106781186547524400844362104849d0 rkB(1) = -.707106781186547524400844362104849d0 rkB(2) = 1.707106781186547524400844362104849d0 rkBhat(1)= .6666666666666666666666666666666667d0 rkBhat(2)= .3333333333333333333333333333333333d0 rkC(1) = 1.707106781186547524400844362104849d0 rkC(2) = 1.0d0 ! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} rkD(1) = 0.0d0 rkD(2) = 1.0d0 ! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} rkE(1) = -.4714045207910316829338962414032326d0 rkE(2) = .8047378541243650162672295747365659d0 ! Local order of Err estimate rkElo = 2 ! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} rkTheta(2,1) = -.414213562373095048801688724209698d0 ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = .5857864376269049511983112757903019d0 END SUBROUTINE Sdirk2b !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Sdirk3a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ sdMethod = S3A ! Number of stages rkS = 3 ! Method coefficients rkGamma = .2113248654051871177454256097490213d0 rkA(1,1) = .2113248654051871177454256097490213d0 rkA(2,1) = .2113248654051871177454256097490213d0 rkA(2,2) = .2113248654051871177454256097490213d0 rkA(3,1) = .2113248654051871177454256097490213d0 rkA(3,2) = .5773502691896257645091487805019573d0 rkA(3,3) = .2113248654051871177454256097490213d0 rkB(1) = .2113248654051871177454256097490213d0 rkB(2) = .5773502691896257645091487805019573d0 rkB(3) = .2113248654051871177454256097490213d0 rkBhat(1)= .2113248654051871177454256097490213d0 rkBhat(2)= .6477918909913548037576239837516312d0 rkBhat(3)= .1408832436034580784969504064993475d0 rkC(1) = .2113248654051871177454256097490213d0 rkC(2) = .4226497308103742354908512194980427d0 rkC(3) = 1.d0 ! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i} rkD(1) = 0.d0 rkD(2) = 0.d0 rkD(3) = 1.d0 ! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i} rkE(1) = 0.9106836025229590978424821138352906d0 rkE(2) = -1.244016935856292431175815447168624d0 rkE(3) = 0.3333333333333333333333333333333333d0 ! Local order of Err estimate rkElo = 2 ! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j} rkTheta(2,1) = 1.0d0 rkTheta(3,1) = -1.732050807568877293527446341505872d0 rkTheta(3,2) = 2.732050807568877293527446341505872d0 ! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j} rkAlpha(2,1) = 2.0d0 rkAlpha(3,1) = -12.92820323027550917410978536602349d0 rkAlpha(3,2) = 8.83012701892219323381861585376468d0 END SUBROUTINE Sdirk3a !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ END SUBROUTINE SdirkTLM ! AND ALL ITS INTERNAL PROCEDURES !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE FUN_CHEM( T, Y, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE KPP_ROOT_Parameters, ONLY: NVAR USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST USE KPP_ROOT_Function, ONLY: Fun USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO KPP_REAL :: T, Told KPP_REAL :: Y(NVAR), P(NVAR) Told = TIME TIME = T CALL Update_SUN() CALL Update_RCONST() CALL Fun( Y, FIX, RCONST, P ) TIME = Told END SUBROUTINE FUN_CHEM !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE JAC_CHEM( T, Y, JV ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST USE KPP_ROOT_Jacobian, ONLY: Jac_SP,LU_IROW,LU_ICOL USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO KPP_REAL :: T, Told KPP_REAL :: Y(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR) INTEGER :: i, j #else KPP_REAL :: JV(LU_NONZERO) #endif Told = TIME TIME = T CALL Update_SUN() CALL Update_RCONST() #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JS) DO j=1,NVAR DO i=1,NVAR JV(i,j) = 0.0D0 END DO END DO DO i=1,LU_NONZERO JV(LU_IROW(i),LU_ICOL(i)) = JS(i) END DO #else CALL Jac_SP(Y, FIX, RCONST, JV) #endif TIME = Told END SUBROUTINE JAC_CHEM END MODULE KPP_ROOT_Integrator