! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! Numerical Integrator (Time-Stepping) File ! ! Generated by KPP-2.2 symbolic chemistry Kinetics PreProcessor ! (http://www.cs.vt.edu/~asandu/Software/KPP) ! KPP is distributed under GPL, the general public licence ! (http://www.gnu.org/copyleft/gpl.html) ! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa ! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech ! With important contributions from: ! M. Damian, Villanova University, USA ! R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany ! ! File : gckpp_adj_Integrator.f90 ! Time : Tue May 14 19:43:54 2013 ! Working directory : /home/daven/kpp-2.2.1/GC_KPP ! Equation file : gckpp_adj.kpp ! Output root filename : gckpp_adj ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! INTEGRATE - Integrator routine ! Arguments : ! TIN - Start Time for Integration ! TOUT - End Time for Integration ! ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! Discrete adjoints of Rosenbrock, ! ! for several Rosenbrock methods: ! ! * Ros2 ! ! * Ros3 ! ! * Ros4 ! ! * Rodas3 ! ! * Rodas4 ! ! By default the code employs the KPP sparse linear algebra routines ! ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! ! ! ! (C) Adrian Sandu, August 2004 ! ! Virginia Polytechnic Institute and State University ! ! Contact: sandu@cs.vt.edu ! ! Revised by Philipp Miehe and Adrian Sandu, May 2006 ! ! Revised by Adrian Sandu, March 2008: ! ! added sensitivity w.r.t. rate coefficients, following D.K. Henze ! ! ! This implementation is part of KPP - the Kinetic PreProcessor ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! MODULE gckpp_adj_Integrator USE gckpp_adj_Precision USE gckpp_adj_Parameters USE gckpp_adj_Global USE gckpp_adj_LinearAlgebra USE gckpp_adj_Rates USE gckpp_adj_Function USE gckpp_adj_Jacobian USE gckpp_adj_Hessian USE gckpp_adj_Util IMPLICIT NONE PUBLIC SAVE !~~~> Statistics on the work performed by the Rosenbrock 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, & Nierr=20 CONTAINS ! Routines in the module gckpp_adj_Integrator !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, Lambda_R, & TIN, TOUT, ATOL_adj, RTOL_adj, & ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Y - Concentrations REAL(kind=dp) :: Y(NVAR) !~~~> NADJ - No. of cost functionals for which adjoints ! are evaluated simultaneously ! If single cost functional is considered (like in ! most applications) simply set NADJ = 1 INTEGER NADJ !~~~> Lambda - Sensitivities w.r.t. concentrations ! Note: Lambda (1:NVAR,j) contains sensitivities of ! the j-th cost functional w.r.t. Y(1:NVAR), j=1...NADJ REAL(kind=dp) :: Lambda(NVAR,NADJ) !~~~> Lambda_R - Sensitivities w.r.t. rate coefficients ! of reactions JCOEFF(1) ... JCOEFF(NCOEFF) REAL(kind=dp) :: Lambda_R(NCOEFF,NADJ) !~~~> Time interval REAL(kind=dp), INTENT(IN) :: TIN ! TIN - Start Time REAL(kind=dp), INTENT(IN) :: TOUT ! TOUT - End Time !~~~> Tolerances for adjoint calculations ! (used only for full continuous adjoint) REAL(kind=dp), INTENT(IN) :: ATOL_adj(NVAR,NADJ), RTOL_adj(NVAR,NADJ) !~~~> Optional input parameters and statistics INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) REAL(kind=dp), INTENT(IN), OPTIONAL :: RCNTRL_U(20) INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) REAL(kind=dp), INTENT(OUT), OPTIONAL :: RSTATUS_U(20) REAL(kind=dp) :: RCNTRL(20), RSTATUS(20) INTEGER :: ICNTRL(20), ISTATUS(20), IERR INTEGER, SAVE :: Ntotal ICNTRL(1:20) = 0 RCNTRL(1:20) = 0.0_dp ISTATUS(1:20) = 0 RSTATUS(1:20) = 0.0_dp !~~~> fine-tune the integrator: ! ICNTRL(1) = 0 ! 0 = non-autonomous, 1 = autonomous ! ICNTRL(2) = 1 ! 0 = scalar, 1 = vector tolerances ! RCNTRL(3) = STEPMIN ! starting step ! ICNTRL(3) = 5 ! choice of the method for forward integration ! ICNTRL(6) = 1 ! choice of the method for continuous adjoint ! ICNTRL(7) = 2 ! 1=none, 2=discrete, 3=full continuous, 4=simplified continuous adjoint ! ICNTRL(8) = 1 ! Save fwd LU factorization: 0 = *don't* save, 1 = save ! if optional parameters are given, and if they are >=0, then they 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 CALL RosenbrockADJ(Y, NADJ, Lambda, Lambda_R, & TIN, TOUT, & ATOL, RTOL, ATOL_adj, RTOL_adj, & RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR) !~~~> Debug option: show number of steps ! Ntotal = Ntotal + ISTATUS(Nstp) ! WRITE(6,777) ISTATUS(Nstp),Ntotal,VAR(ind_O3),VAR(ind_NO2) !777 FORMAT('NSTEPS=',I6,' (',I6,') O3=',E24.14,' NO2=',E24.14) IF (IERR < 0) THEN print *,'RosenbrockADJ: Unsucessful step at T=', & TIN,' (IERR=',IERR,')' ! Now record IERR value in RSATUS(Niere) (dkh, 07/05/06) ISTATUS(Nierr) = IERR END IF STEPMIN = RSTATUS(Nhexit) ! if optional parameters are given for output ! copy to them to return information IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:) IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:) END SUBROUTINE INTEGRATE_ADJ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, Lambda_R, & Tstart, Tend, & AbsTol, RelTol, AbsTol_adj, RelTol_adj, & RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! ! ADJ = Adjoint of the Tangent Linear Model of a Rosenbrock Method ! ! Solves the system y'=F(t,y) using a RosenbrockADJ method defined by: ! ! G = 1/(H*gamma(1)) - Jac(t0,Y0) ! T_i = t0 + Alpha(i)*H ! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j ! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j + ! gamma(i)*dF/dT(t0, Y0) ! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j ! ! For details on RosenbrockADJ methods and their implementation consult: ! E. Hairer and G. Wanner ! "Solving ODEs II. Stiff and differential-algebraic problems". ! Springer series in computational mathematics, Springer-Verlag, 1996. ! The codes contained in the book inspired this implementation. ! ! (C) Adrian Sandu, August 2004 ! 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=Tstart) ! NADJ -> dimension of linearized system, ! i.e. the number of sensitivity coefficients !- Lambda(NVAR,NADJ) -> vector of initial sensitivity conditions (at T=Tstart) !- [Tstart,Tend] = time range of integration ! (if Tstart>Tend the integration is performed backwards in time) !- RelTol, AbsTol = user precribed accuracy !- SUBROUTINE Fun( T, Y, Ydot ) = ODE function, ! returns Ydot = Y' = F(T,Y) !- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function, ! returns Jcb = dF/dY !- ICNTRL(1:10) = integer inputs parameters !- RCNTRL(1:10) = real inputs parameters !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> OUTPUT ARGUMENTS: ! !- Y(NVAR) -> vector of final states (at T->Tend) !- Lambda(NVAR,NADJ) -> vector of final sensitivities (at T=Tend) !- ICNTRL(11:20) -> integer output parameters !- RCNTRL(11:20) -> real output parameters !- IERR -> job status upon return ! - succes (positive value) or failure (negative value) - ! = 1 : Success ! = -1 : Improper value for maximal no of steps ! = -2 : Selected RosenbrockADJ method not implemented ! = -3 : Hmin/Hmax/Hstart must be positive ! = -4 : FacMin/FacMax/FacRej must be positive ! = -5 : Improper tolerance values ! = -6 : No of steps exceeds maximum bound ! = -7 : Step size too small ! = -8 : Matrix is repeatedly singular !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> INPUT PARAMETERS: ! ! Note: For input parameters equal to zero the default values of the ! corresponding variables are used. ! ! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS) ! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS) ! ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors ! = 1: AbsTol, RelTol are scalars ! ! ICNTRL(3) -> selection of a particular Rosenbrock method ! = 0 : default method is Rodas3 ! = 1 : method is Ros2 ! = 2 : method is Ros3 ! = 3 : method is Ros4 ! = 4 : method is Rodas3 ! = 5: method is Rodas4 ! ! ICNTRL(4) -> maximum number of integration steps ! For ICNTRL(4)=0) the default value of BUFSIZE is used ! ! ICNTRL(6) -> selection of a particular Rosenbrock method for the ! continuous adjoint integration - for cts adjoint it ! can be different than the forward method ICNTRL(3) ! Note 1: to avoid interpolation errors (which can be huge!) ! it is recommended to use only ICNTRL(7) = 2 or 4 ! Note 2: the performance of the full continuous adjoint ! strongly depends on the forward solution accuracy Abs/RelTol ! ! ICNTRL(7) -> Type of adjoint algorithm ! = 0 : default is discrete adjoint ( of method ICNTRL(3) ) ! plus sensitivity w.r.t. reaction coefficients ! = 1 : no adjoint ! = 2 : discrete adjoint ( of method ICNTRL(3) ) ! plus sensitivity w.r.t. reaction coefficients ! = 3 : fully adaptive continuous adjoint ( with method ICNTRL(6) ) ! = 4 : simplified continuous adjoint ( with method ICNTRL(6) ) ! ! ICNTRL(8) -> checkpointing the LU factorization at each step: ! ICNTRL(8)=0 : do *not* save LU factorization (the default) ! ICNTRL(8)=1 : save LU factorization ! Note: if ICNTRL(7)=1 the LU factorization is *not* saved ! !~~~> Real input 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) ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! !~~~> OUTPUT PARAMETERS: ! ! Note: each call to RosenbrockADJ adds the corrent 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 exit ! For multiple restarts, use Hexit as Hstart in the following run !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Arguments REAL(kind=dp), INTENT(INOUT) :: Y(NVAR) INTEGER, INTENT(IN) :: NADJ REAL(kind=dp), INTENT(INOUT) :: Lambda(NVAR,NADJ) REAL(kind=dp), INTENT(INOUT) :: Lambda_R(NCOEFF,NADJ) REAL(kind=dp), INTENT(IN) :: Tstart,Tend REAL(kind=dp), INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR) REAL(kind=dp), INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ) INTEGER, INTENT(IN) :: ICNTRL(20) REAL(kind=dp), INTENT(IN) :: RCNTRL(20) INTEGER, INTENT(INOUT) :: ISTATUS(20) REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20) INTEGER, INTENT(OUT) :: IERR !~~~> Parameters of the Rosenbrock method, up to 6 stages INTEGER :: ros_S, rosMethod INTEGER, PARAMETER :: RS2=1, RS3=2, RS4=3, RD3=4, RD4=5 REAL(kind=dp) :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), & ros_Alpha(6), ros_Gamma(6), ros_ELO LOGICAL :: ros_NewF(6) CHARACTER(LEN=12) :: ros_Name !~~~> Types of Adjoints Implemented INTEGER, PARAMETER :: Adj_none = 1, Adj_discrete = 2, & Adj_continuous = 3, Adj_simple_continuous = 4 !~~~> Checkpoints in memory ! Can make this much smaller (dkh, 01/06/10) !INTEGER, PARAMETER :: bufsize = 200000 INTEGER, PARAMETER :: bufsize = 4000 ! Need to make stack_ptr THREADPRIVATE (dkh, 07/31/09) !INTEGER :: stack_ptr = 0 ! last written entry REAL(kind=dp), DIMENSION(:), POINTER :: chk_H, chk_T REAL(kind=dp), DIMENSION(:,:), POINTER :: chk_Y, chk_K, chk_J REAL(kind=dp), DIMENSION(:,:), POINTER :: chk_dY, chk_d2Y !~~~> Local variables REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe REAL(kind=dp) :: Hmin, Hmax, Hstart REAL(kind=dp) :: Texit INTEGER :: i, UplimTol, Max_no_steps INTEGER :: IERR_SAVE INTEGER :: AdjointType, CadjMethod LOGICAL :: Autonomous, VectorTol, SaveLU !~~~> Parameters REAL(kind=dp), PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 REAL(kind=dp), PARAMETER :: DeltaMin = 1.0d-5 !~~~> Initialize statistics ISTATUS(1:20) = 0 RSTATUS(1:20) = ZERO !~~~> Autonomous or time dependent ODE. Default is time dependent. Autonomous = .NOT.(ICNTRL(1) == 0) !~~~> 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 VectorTol = .TRUE. UplimTol = NVAR ELSE VectorTol = .FALSE. UplimTol = 1 END IF !~~~> Initialize the particular Rosenbrock method selected SELECT CASE (ICNTRL(3)) CASE (1) CALL Ros2 CASE (2) CALL Ros3 CASE (3) CALL Ros4 CASE (0,4) CALL Rodas3 CASE (5) CALL Rodas4 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) RETURN END SELECT !~~~> The maximum number of steps admitted IF (ICNTRL(4) == 0) THEN Max_no_steps = bufsize - 1 ELSEIF (Max_no_steps > 0) THEN Max_no_steps=ICNTRL(4) ELSE PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4) CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR) RETURN END IF !~~~> The particular Rosenbrock method chosen for integrating the cts adjoint IF (ICNTRL(6) == 0) THEN CadjMethod = 4 ELSEIF ( (ICNTRL(6) >= 1).AND.(ICNTRL(6) <= 5) ) THEN CadjMethod = ICNTRL(6) ELSE PRINT * , 'Unknown CADJ Rosenbrock method: ICNTRL(6)=', CadjMethod CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) RETURN END IF !~~~> Discrete or continuous adjoint formulation IF ( ICNTRL(7) == 0 ) THEN AdjointType = Adj_discrete ELSEIF ( (ICNTRL(7) >= 1).AND.(ICNTRL(7) <= 4) ) THEN AdjointType = ICNTRL(7) ELSE PRINT * , 'User-selected adjoint type: ICNTRL(7)=', AdjointType CALL ros_ErrorMsg(-9,Tstart,ZERO,IERR) RETURN END IF !~~~> Save or not the forward LU factorization SaveLU = (ICNTRL(8) /= 0) !~~~> 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 Hmin: RCNTRL(1)=', RCNTRL(1) CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR) RETURN END IF !~~~> Upper bound on the step size: (positive value) IF (RCNTRL(2) == ZERO) THEN Hmax = ABS(Tend-Tstart) ELSEIF (RCNTRL(2) > ZERO) THEN Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart)) ELSE PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2) CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR) RETURN END IF !~~~> Starting step size: (positive value) IF (RCNTRL(3) == ZERO) THEN Hstart = MAX(Hmin,DeltaMin) ELSEIF (RCNTRL(3) > ZERO) THEN Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart)) ELSE PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR) RETURN END IF !~~~> Step size can be changed s.t. FacMin < Hnew/Hold < FacMax IF (RCNTRL(4) == ZERO) THEN FacMin = 0.2d0 ELSEIF (RCNTRL(4) > ZERO) THEN FacMin = RCNTRL(4) ELSE PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4) CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR) RETURN END IF IF (RCNTRL(5) == ZERO) THEN FacMax = 6.0d0 ELSEIF (RCNTRL(5) > ZERO) THEN FacMax = RCNTRL(5) ELSE PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5) CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR) RETURN END IF !~~~> FacRej: Factor to decrease step after 2 succesive rejections IF (RCNTRL(6) == ZERO) THEN FacRej = 0.1d0 ELSEIF (RCNTRL(6) > ZERO) THEN FacRej = RCNTRL(6) ELSE PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6) CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR) RETURN END IF !~~~> FacSafe: Safety Factor in the computation of new step size IF (RCNTRL(7) == ZERO) THEN FacSafe = 0.9d0 ELSEIF (RCNTRL(7) > ZERO) THEN FacSafe = RCNTRL(7) ELSE PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7) CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR) RETURN END IF !~~~> Check if tolerances are reasonable DO i=1,UplimTol IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.d0*Roundoff) & .OR. (RelTol(i) >= 1.0d0) ) THEN PRINT * , ' AbsTol(',i,') = ',AbsTol(i) PRINT * , ' RelTol(',i,') = ',RelTol(i) CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR) RETURN END IF END DO !~~~> Allocate checkpoint space or open checkpoint files IF (AdjointType == Adj_discrete) THEN CALL ros_AllocateDBuffers( ros_S ) ELSEIF ( (AdjointType == Adj_continuous).OR. & (AdjointType == Adj_simple_continuous) ) THEN CALL ros_AllocateCBuffers END IF !~~~> CALL Forward Rosenbrock method CALL ros_FwdInt(Y,Tstart,Tend,Texit, & AbsTol, RelTol, & ! Error indicator IERR) !!$ PRINT*,'FORWARD STATISTICS' !!$ PRINT*,'Step=',Nstp,' Acc=',Nacc, & !!$ ' Rej=',Nrej, ' Singular=',Nsng ! Now contiue to avoid having adjoint arrays become corrupt (dkh, 07/08/11, adj32_004) !~~~> If Forward integration failed return ! IF (IERR<0) RETURN ! but save a copy of the error stat IERR_SAVE = IERR !~~~> Initialize the particular Rosenbrock method for continuous adjoint IF ( (AdjointType == Adj_continuous).OR. & (AdjointType == Adj_simple_continuous) ) THEN SELECT CASE (CadjMethod) CASE (1) CALL Ros2 CASE (2) CALL Ros3 CASE (3) CALL Ros4 CASE (4) CALL Rodas3 CASE (5) CALL Rodas4 CASE DEFAULT PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', ICNTRL(3) CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR) RETURN END SELECT END IF SELECT CASE (AdjointType) CASE (Adj_discrete) CALL ros_DadjRateInt ( & NADJ, Lambda, Lambda_R, & Tstart, Tend, Texit, & IERR ) CASE (Adj_continuous) CALL ros_CadjInt ( & NADJ, Lambda, & Tend, Tstart, Texit, & AbsTol_adj, RelTol_adj, & IERR ) CASE (Adj_simple_continuous) CALL ros_SimpleCadjInt ( & NADJ, Lambda, & Tstart, Tend, Texit, & IERR ) END SELECT ! AdjointType !!$ PRINT*,'ADJOINT STATISTICS' !!$ PRINT*,'Step=',Nstp,' Acc=',Nacc, & !!$ ' Rej=',Nrej, ' Singular=',Nsng !~~~> Free checkpoint space or close checkpoint files IF (AdjointType == Adj_discrete) THEN CALL ros_FreeDBuffers ELSEIF ( (AdjointType == Adj_continuous) .OR. & (AdjointType == Adj_simple_continuous) ) THEN CALL ros_FreeCBuffers END IF ! replace with original error stat from fwd int dkh IERR = IERR_SAVE !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CONTAINS ! Procedures internal to RosenbrockADJ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_AllocateDBuffers( S ) !~~~> Allocate buffer space for discrete adjoint !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i, S ALLOCATE( chk_H(bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer H'; STOP END IF ALLOCATE( chk_T(bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer T'; STOP END IF ALLOCATE( chk_Y(NVAR*S,bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer Y'; STOP END IF ALLOCATE( chk_K(NVAR*S,bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer K'; STOP END IF IF (SaveLU) THEN #ifdef FULL_ALGEBRA ALLOCATE( chk_J(NVAR*NVAR,bufsize), STAT=i ) #else ALLOCATE( chk_J(LU_NONZERO,bufsize), STAT=i ) #endif IF (i/=0) THEN PRINT*,'Failed allocation of buffer J'; STOP END IF END IF END SUBROUTINE ros_AllocateDBuffers !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_FreeDBuffers !~~~> Dallocate buffer space for discrete adjoint !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i DEALLOCATE( chk_H, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer H'; STOP END IF DEALLOCATE( chk_T, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer T'; STOP END IF DEALLOCATE( chk_Y, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer Y'; STOP END IF DEALLOCATE( chk_K, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer K'; STOP END IF IF (SaveLU) THEN DEALLOCATE( chk_J, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer J'; STOP END IF END IF END SUBROUTINE ros_FreeDBuffers !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_AllocateCBuffers !~~~> Allocate buffer space for continuous adjoint !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i ALLOCATE( chk_H(bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer H'; STOP END IF ALLOCATE( chk_T(bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer T'; STOP END IF ALLOCATE( chk_Y(NVAR,bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer Y'; STOP END IF ALLOCATE( chk_dY(NVAR,bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer dY'; STOP END IF ALLOCATE( chk_d2Y(NVAR,bufsize), STAT=i ) IF (i/=0) THEN PRINT*,'Failed allocation of buffer d2Y'; STOP END IF END SUBROUTINE ros_AllocateCBuffers !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_FreeCBuffers !~~~> Dallocate buffer space for continuous adjoint !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: i DEALLOCATE( chk_H, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer H'; STOP END IF DEALLOCATE( chk_T, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer T'; STOP END IF DEALLOCATE( chk_Y, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer Y'; STOP END IF DEALLOCATE( chk_dY, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer dY'; STOP END IF DEALLOCATE( chk_d2Y, STAT=i ) IF (i/=0) THEN PRINT*,'Failed deallocation of buffer d2Y'; STOP END IF END SUBROUTINE ros_FreeCBuffers !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_DPush( S, T, H, Ystage, K, E, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Saves the next trajectory snapshot for discrete adjoints !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: S ! no of stages REAL(kind=dp) :: T, H, Ystage(NVAR*S), K(NVAR*S) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA REAL(kind=dp) :: E(NVAR,NVAR) #else REAL(kind=dp) :: E(LU_NONZERO) #endif stack_ptr = stack_ptr + 1 IF ( stack_ptr > bufsize ) THEN PRINT*,'Push failed: buffer overflow' STOP END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T !CALL WCOPY(NVAR*S,Ystage,1,chk_Y(1,stack_ptr),1) !CALL WCOPY(NVAR*S,K,1,chk_K(1,stack_ptr),1) chk_Y(1:NVAR*S,stack_ptr) = Ystage(1:NVAR*S) chk_K(1:NVAR*S,stack_ptr) = K(1:NVAR*S) IF (SaveLU) THEN #ifdef FULL_ALGEBRA chk_J(1:NVAR,1:NVAR,stack_ptr) = E(1:NVAR,1:NVAR) chk_P(1:NVAR,stack_ptr) = P(1:NVAR) #else chk_J(1:LU_NONZERO,stack_ptr) = E(1:LU_NONZERO) #endif END IF END SUBROUTINE ros_DPush !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_DPop( S, T, H, Ystage, K, E, P ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Retrieves the next trajectory snapshot for discrete adjoints !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ INTEGER :: S ! no of stages REAL(kind=dp) :: T, H, Ystage(NVAR*S), K(NVAR*S) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA REAL(kind=dp) :: E(NVAR,NVAR) #else REAL(kind=dp) :: E(LU_NONZERO) #endif IF ( stack_ptr <= 0 ) THEN PRINT*,'Pop failed: empty buffer' STOP END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) !CALL WCOPY(NVAR*S,chk_Y(1,stack_ptr),1,Ystage,1) !CALL WCOPY(NVAR*S,chk_K(1,stack_ptr),1,K,1) Ystage(1:NVAR*S) = chk_Y(1:NVAR*S,stack_ptr) K(1:NVAR*S) = chk_K(1:NVAR*S,stack_ptr) !CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1) IF (SaveLU) THEN #ifdef FULL_ALGEBRA E(1:NVAR,1:NVAR) = chk_J(1:NVAR,1:NVAR,stack_ptr) P(1:NVAR) = chk_P(1:NVAR,stack_ptr) #else E(1:LU_NONZERO) = chk_J(1:LU_NONZERO,stack_ptr) #endif END IF stack_ptr = stack_ptr - 1 END SUBROUTINE ros_DPop !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_CPush( T, H, Y, dY, d2Y ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Saves the next trajectory snapshot for discrete adjoints !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL(kind=dp) :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR) stack_ptr = stack_ptr + 1 IF ( stack_ptr > bufsize ) THEN PRINT*,'Push failed: buffer overflow' STOP END IF chk_H( stack_ptr ) = H chk_T( stack_ptr ) = T !CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1) !CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1) !CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1) chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR) chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR) chk_d2Y(1:NVAR,stack_ptr) = d2Y(1:NVAR) END SUBROUTINE ros_CPush !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_CPop( T, H, Y, dY, d2Y ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Retrieves the next trajectory snapshot for discrete adjoints !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL(kind=dp) :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR) IF ( stack_ptr <= 0 ) THEN PRINT*,'Pop failed: empty buffer' STOP END IF H = chk_H( stack_ptr ) T = chk_T( stack_ptr ) !CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1) !CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1) !CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1) Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr) dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr) d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr) stack_ptr = stack_ptr - 1 END SUBROUTINE ros_CPop !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_ErrorMsg(Code,T,H,IERR) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Handles all error messages !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL(kind=dp), INTENT(IN) :: T, H INTEGER, INTENT(IN) :: Code INTEGER, INTENT(OUT) :: IERR IERR = Code PRINT * , & 'Forced exit from RosenbrockADJ due to the following error:' SELECT CASE (Code) CASE (-1) PRINT * , '--> Improper value for maximal no of steps' CASE (-2) PRINT * , '--> Selected RosenbrockADJ method not implemented' 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 buffer bound' CASE (-7) PRINT * , '--> Step size too small: T + 10*H = T', & ' or H < Roundoff' CASE (-8) PRINT * , '--> Matrix is repeatedly singular' CASE (-9) PRINT * , '--> Improper type of adjoint selected' CASE DEFAULT PRINT *, 'Unknown Error code: ', Code END SELECT PRINT *, 'T=', T, 'and H=', H END SUBROUTINE ros_ErrorMsg !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_FwdInt (Y, & Tstart, Tend, T, & AbsTol, RelTol, & !~~~> Error indicator IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the implementation of a generic RosenbrockADJ method ! defined by ros_S (no of stages) ! and its coefficients ros_{A,C,M,E,Alpha,Gamma} !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input: the initial condition at Tstart; Output: the solution at T REAL(kind=dp), INTENT(INOUT) :: Y(NVAR) !~~~> Input: integration interval REAL(kind=dp), INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) REAL(kind=dp), INTENT(OUT) :: T !~~~> Input: tolerances REAL(kind=dp), INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR) !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables REAL(kind=dp) :: Ynew(NVAR), Fcn0(NVAR), Fcn(NVAR) REAL(kind=dp) :: K(NVAR*ros_S), dFdT(NVAR) REAL(kind=dp), DIMENSION(:), POINTER :: Ystage #ifdef FULL_ALGEBRA REAL(kind=dp) :: Jac0(NVAR,NVAR), Ghimj(NVAR,NVAR) #else REAL(kind=dp) :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO) #endif REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau REAL(kind=dp) :: Err, Yerr(NVAR) INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage LOGICAL :: RejectLastH, RejectMoreH, Singular !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Allocate stage vector buffer if needed IF (AdjointType == Adj_discrete) THEN ALLOCATE(Ystage(NVAR*ros_S), STAT=i) ! Uninitialized Ystage may lead to NaN on some compilers Ystage = 0.0d0 IF (i/=0) THEN PRINT*,'Allocation of Ystage failed' STOP END IF END IF !~~~> Initial preparations T = Tstart RSTATUS(Nhexit) = ZERO H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) ) IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin IF (Tend >= Tstart) THEN Direction = +1 ELSE Direction = -1 END IF H = Direction*H RejectLastH=.FALSE. RejectMoreH=.FALSE. !~~~> Time loop begins below TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff*abs(Tend) <= ZERO) & .OR. (Direction < 0).AND.((Tend-T)+Roundoff*abs(Tend) <= ZERO) ) ! Added *abs(Tend) by KS, A.Sandu for boundary cases IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps CALL ros_ErrorMsg(-6,T,H,IERR) RETURN END IF IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small CALL ros_ErrorMsg(-7,T,H,IERR) RETURN END IF !~~~> Limit H if necessary to avoid going beyond Tend RSTATUS(Nhexit) = H H = MIN(H,ABS(Tend-T)) !~~~> Compute the function at current time CALL FunTemplate(T,Y,Fcn0) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 !~~~> Compute the function derivative with respect to T IF (.NOT.Autonomous) THEN CALL ros_FunTimeDerivative ( T, Roundoff, Y, & Fcn0, dFdT ) END IF !~~~> Compute the Jacobian at current time CALL JacTemplate(T,Y,Jac0) ISTATUS(Njac) = ISTATUS(Njac) + 1 !~~~> Repeat step calculation until current step accepted UntilAccepted: DO CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), & Jac0,Ghimj,Pivot,Singular) IF (Singular) THEN ! More than 5 consecutive failed decompositions CALL ros_ErrorMsg(-8,T,H,IERR) RETURN END IF !~~~> Compute the stages Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR) ioffset = NVAR*(istage-1) ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN CALL WCOPY(NVAR,Fcn0,1,Fcn,1) IF (AdjointType == Adj_discrete) THEN ! Save stage solution ! CALL WCOPY(NVAR,Y,1,Ystage(1),1) Ystage(1:NVAR) = Y(1:NVAR) CALL WCOPY(NVAR,Y,1,Ynew,1) END IF ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN CALL WCOPY(NVAR,Y,1,Ynew,1) DO j = 1, istage-1 CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & K(NVAR*(j-1)+1),1,Ynew,1) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL FunTemplate(Tau,Ynew,Fcn) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 END IF ! if istage == 1 elseif ros_NewF(istage) ! save stage solution every time even if ynew is not updated IF ( ( istage > 1 ).AND.(AdjointType == Adj_discrete) ) THEN ! CALL WCOPY(NVAR,Ynew,1,Ystage(ioffset+1),1) Ystage(ioffset+1:ioffset+NVAR) = Ynew(1:NVAR) END IF CALL WCOPY(NVAR,Fcn,1,K(ioffset+1),1) DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1),1,K(ioffset+1),1) END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1),1) END IF CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1)) END DO Stage !~~~> Compute the new solution CALL WCOPY(NVAR,Y,1,Ynew,1) DO j=1,ros_S CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1),1,Ynew,1) END DO !~~~> Compute the error estimation CALL WSCAL(NVAR,ZERO,Yerr,1) DO j=1,ros_S CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1),1,Yerr,1) END DO Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol ) !~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO))) Hnew = H*Fac !~~~> Check the error magnitude and adjust step size ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 IF (AdjointType == Adj_discrete) THEN ! Save current state CALL ros_DPush( ros_S, T, H, Ystage, K, Ghimj, Pivot ) ELSEIF ( (AdjointType == Adj_continuous) .OR. & (AdjointType == Adj_simple_continuous) ) THEN #ifdef FULL_ALGEBRA K = MATMUL(Jac0,Fcn0) #else CALL Jac_SP_Vec( Jac0, Fcn0, K(1) ) #endif IF (.NOT. Autonomous) THEN CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1) ) END IF CALL WCOPY(NVAR,Ynew,1,Y,1) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step Hnew = MIN(Hnew,H) END IF RSTATUS(Nhexit) = H RSTATUS(Nhnew) = Hnew RSTATUS(Ntexit) = T RejectLastH = .FALSE. RejectMoreH = .FALSE. H = Hnew EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED ELSE !~~~> Reject step IF (RejectMoreH) THEN Hnew = H*FacRej END IF RejectMoreH = RejectLastH RejectLastH = .TRUE. H = Hnew IF (ISTATUS(Nacc) >= 1) THEN ISTATUS(Nrej) = ISTATUS(Nrej) + 1 END IF END IF ! Err <= 1 END DO UntilAccepted END DO TimeLoop !~~~> Save last state: only needed for continuous adjoint IF ( (AdjointType == Adj_continuous) .OR. & (AdjointType == Adj_simple_continuous) ) THEN CALL FunTemplate(T,Y,Fcn0) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 CALL JacTemplate(T,Y,Jac0) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA K = MATMUL(Jac0,Fcn0) #else CALL Jac_SP_Vec( Jac0, Fcn0, K(1) ) #endif IF (.NOT. Autonomous) THEN CALL ros_FunTimeDerivative ( T, Roundoff, Y, & Fcn0, dFdT ) CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1) ) !~~~> Deallocate stage buffer: only needed for discrete adjoint ELSEIF (AdjointType == Adj_discrete) THEN DEALLOCATE(Ystage, STAT=i) IF (i/=0) THEN PRINT*,'Deallocation of Ystage failed' STOP END IF END IF !~~~> Succesful exit IERR = 1 !~~~> The integration was successful END SUBROUTINE ros_FwdInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_DadjRateInt ( & NADJ, Lambda, Lambda_R, & Tstart, Tend, T, & !~~~> Error indicator IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the implementation of a generic RosenbrockADJ method ! defined by ros_S (no of stages) ! and its coefficients ros_{A,C,M,E,Alpha,Gamma} ! The adjoint sensitivity of the solution with respect to NCOEFF selected ! reaction rate coefficients JCOEFF(1:NCOEFF) is also included ! Note: works only for autonomous systems, with fixed (in time) RCOEFF ! Based on the implementation of Daven K. Henze !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Contains dFun_dRcoeff, dJac_dRcoeff USE gckpp_adj_STOICHIOM ! added logical switch (tww, 05/08/12) USE LOGICAL_ADJ_MOD, ONLY : LADJ_RRATE IMPLICIT NONE !~~~> Input: the initial condition at Tstart; Output: the solution at T INTEGER, INTENT(IN) :: NADJ !~~~> First order adjoint REAL(kind=dp), INTENT(INOUT) :: Lambda(NVAR,NADJ) !!~~~> Input: integration interval REAL(kind=dp), INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) REAL(kind=dp), INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables REAL(kind=dp) :: Ystage(NVAR*ros_S), K(NVAR*ros_S) REAL(kind=dp) :: U(NVAR*ros_S,NADJ), V(NVAR*ros_S,NADJ) #ifdef FULL_ALGEBRA REAL(kind=dp), DIMENSION(NVAR,NVAR) :: Jac, dJdT, Ghimj #else REAL(kind=dp), DIMENSION(LU_NONZERO) :: Jac, dJdT, Ghimj #endif REAL(kind=dp) :: Hes0(NHESS) REAL(kind=dp) :: Tmp(NVAR), Tmp2(NVAR) REAL(kind=dp) :: H, HC, HA, Tau INTEGER :: Pivot(NVAR), Direction INTEGER :: i, j, m, istage, istart, jstart !~~~> Local parameters REAL(kind=dp), PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 REAL(kind=dp), PARAMETER :: DeltaMin = 1.0d-5 !~~~> Sensitivities w.r.t. reacton coefficients REAL(kind=dp), INTENT(OUT) :: Lambda_R(NCOEFF,NADJ) INTEGER :: icoeff, vstart REAL(kind=dp) :: DFDR(NVAR*NCOEFF) REAL(kind=dp) :: DJDR(NVAR*NCOEFF) REAL(kind=dp) :: DJDR_O3dep(NVAR) REAL(kind=dp) :: V_R(NCOEFF*ros_S,NADJ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (.NOT.Autonomous) THEN PRINT*,'ERROR: ros_DadjRateInt cannot handle NON-AUTONOMOUS systems' STOP END IF IF (Tend >= Tstart) THEN Direction = +1 ELSE Direction = -1 END IF !~~~> dkh Initialize added to 0d0 at the start of each loop Lambda_R(:,:) = ZERO DFDR(:) = ZERO DJDR(:) = ZERO DJDR_O3dep(:) = ZERO V_R(:,:) = ZERO !~~~> Time loop begins below TimeLoop: DO WHILE ( stack_ptr > 0 ) !~~~> Recover checkpoints for stage values and vectors CALL ros_DPop( ros_S, T, H, Ystage, K, Ghimj, Pivot ) ! ISTATUS(Nstp) = ISTATUS(Nstp) + 1 !~~~> Compute LU decomposition IF (.NOT.SaveLU) THEN CALL JacTemplate(T,Ystage(1),Ghimj) ISTATUS(Njac) = ISTATUS(Njac) + 1 Tau = ONE/(Direction*H*ros_Gamma(1)) #ifdef FULL_ALGEBRA Ghimj(1:NVAR,1:NVAR) = -Ghimj(1:NVAR,1:NVAR) DO i=1,NVAR Ghimj(i,i) = Ghimj(i,i)+Tau END DO #else CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau END DO #endif CALL ros_Decomp( Ghimj, Pivot, j ) END IF !~~~> Compute Hessian at the beginning of the interval CALL HessTemplate(T,Ystage(1),Hes0) !~~~> Compute the stages Stage: DO istage = ros_S, 1, -1 !~~~> Current istage first entry istart = NVAR*(istage-1) + 1 !~~~> Compute U DO m = 1,NADJ CALL WCOPY(NVAR,Lambda(1,m),1,U(istart,m),1) CALL WSCAL(NVAR,ros_M(istage),U(istart,m),1) END DO ! m=1:NADJ DO j = istage+1, ros_S jstart = NVAR*(j-1) + 1 HA = ros_A((j-1)*(j-2)/2+istage) HC = ros_C((j-1)*(j-2)/2+istage)/(Direction*H) DO m = 1,NADJ CALL WAXPY(NVAR,HA,V(jstart,m),1,U(istart,m),1) CALL WAXPY(NVAR,HC,U(jstart,m),1,U(istart,m),1) END DO ! m=1:NADJ END DO DO m = 1,NADJ CALL ros_Solve('T', Ghimj, Pivot, U(istart,m)) END DO ! m=1:NADJ !~~~> Compute V Tau = T + ros_Alpha(istage)*Direction*H CALL JacTemplate(Tau,Ystage(istart),Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 ! f_r(Y_istage) CALL dFun_dRcoeff(Ystage(istart:istart+NVAR-1), FIX, NCOEFF, JCOEFF, DFDR) !??CALL dJac_dRcoeff(Ystage(istart:istart+NVAR-1), FIX, K(istart:istart+NVAR-1),NCOEFF, JCOEFF, DJDR) ! don't need for just emissions (dkh, 03/31/10) !! J_r(y_n) x K_istage !CALL dJac_dRcoeff(Ystage(1:NVAR), FIX, K(istart:istart+NVAR-1),NCOEFF, JCOEFF, DJDR) ! added logical switch (tww, 05/08/12) ! Compute for all rxns IF ( LADJ_RRATE ) THEN CALL dJac_dRcoeff(Ystage(1:NVAR), FIX, K(istart:istart+NVAR-1),NCOEFF, JCOEFF, DJDR) ELSE ! Default: Only compute for O3 depostion (needed for shipping emissions), which ! will be the last one. CALL dJac_dRcoeff(Ystage(1:NVAR), FIX, K(istart:istart+NVAR-1),1, JCOEFF(NCOEFF), DJDR_O3dep) ENDIF DO m = 1,NADJ #ifdef FULL_ALGEBRA V(istart:istart+NVAR-1,m) = MATMUL(TRANSPOSE(Jac),U(istart:istart+NVAR-1,m)) #else !---------------------------------------------------------------------- ! Sensitivity w.r.t. reaction rate coefficients (dkh) vstart = NCOEFF*(istage-1) IF ( LADJ_RRATE ) THEN DO icoeff = 1, NCOEFF V_R(vstart+icoeff,m) = ZERO DO j = 1, NVAR !j = DMAP(icoeff) ! += f_r(Y_istage) * U_istage V_R(vstart+icoeff,m) = V_R(vstart+icoeff,m) + DFDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) !! dkh debug !IF ( icoeff == NCOEFF ) THEN ! print*, ' DFDR j = ', DFDR(NVAR*(icoeff-1)+j), j ! print*, ' DJDR j = ', DJDR(NVAR*(icoeff-1)+j), j !ENDIF ! don't need for just emissions (dkh, 03/31/10) !! += ( J_r(y_n) x K_istage )^T * U_istage !V_R(vstart+icoeff,m) = V_R(vstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) ! also comment out calc of DJDR above ! ! added logical switch (tww, 05/08/12) V_R(vstart+icoeff,m) = V_R(vstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) END DO END DO ELSE DO icoeff = 1, NCOEFF-1 V_R(vstart+icoeff,m) = ZERO j = DMAP(icoeff) ! += f_r(Y_istage) * U_istage V_R(vstart+icoeff,m) = V_R(vstart+icoeff,m) + DFDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) ! don't need for just emissions (dkh, 03/31/10) !! += ( J_r(y_n) x K_istage )^T * U_istage !V_R(vstart+icoeff,m) = V_R(vstart+icoeff,m) + DJDR(NVAR*(icoeff-1)+j)*U(istart+j-1,m) END DO ! Now include O3 drydep to account for shipping emissions icoeff = NCOEFF V_R(vstart+icoeff,m) = ZERO V_R(vstart+icoeff,m) = V_R(vstart+icoeff,m) & + DFDR(NVAR*(icoeff-1)+ind_O3 ) * U(istart+ind_O3 -1,m) & + DFDR(NVAR*(icoeff-1)+ind_DRYO3 ) * U(istart+ind_DRYO3 -1,m) & + DFDR(NVAR*(icoeff-1)+ind_DRYDEP) * U(istart+ind_DRYDEP-1,m) & + DJDR_O3dep(ind_O3 ) * U(istart+ind_O3 -1,m) & + DJDR_O3dep(ind_DRYO3 ) * U(istart+ind_DRYO3 -1,m) & + DJDR_O3dep(ind_DRYDEP) * U(istart+ind_DRYDEP-1,m) !! dkh debug !print *, ' DFDR(NVAR*(icoeff-1)+ind_O3 ) = ', DFDR(NVAR*(icoeff-1)+ind_O3 ) !print *, ' DFDR(NVAR*(icoeff-1)+ind_DRYO3 ) = ', DFDR(NVAR*(icoeff-1)+ind_DRYO3 ) !print *, ' DFDR(NVAR*(icoeff-1)+ind_DRYDEP) =', DFDR(NVAR*(icoeff-1)+ind_DRYDEP) !print *, ' DJDR_O3dep(ind_O3 ) = ', DJDR_O3dep(ind_O3 ) !print *, ' DJDR_O3dep(ind_DRYO3 ) = ', DJDR_O3dep(ind_DRYO3 ) !print *, ' DJDR_O3dep(ind_DRYDEP) = ', DJDR_O3dep(ind_DRYDEP) ENDIF !---------------------------------------------------------------------- CALL JacTR_SP_Vec(Jac,U(istart,m),V(istart,m)) #endif END DO ! m=1:NADJ END DO Stage !$$ IF (.NOT.Autonomous) THEN !$$!~~~> Compute the Jacobian derivative with respect to T. !$$! Last "Jac" computed for stage 1 !$$ CALL ros_JacTimeDerivative ( T, Roundoff, Ystage(1), Jac, dJdT ) !$$ END IF !~~~> Compute the new solution !~~~> Compute Lambda_R (dkh) DO istage=1,ros_S vstart = NCOEFF*(istage-1) + 1 DO m = 1,NADJ ! Sum_over_i f_r(Y_i)^T * U_i + (J_r x K_i)^T * U_i CALL WAXPY(NCOEFF,ONE,V_R(vstart,m),1,Lambda_R(1,m),1) END DO ! m=1:NADJ END DO !~~~> Compute Lambda DO istage=1,ros_S istart = NVAR*(istage-1) + 1 DO m = 1,NADJ ! Add V_i CALL WAXPY(NVAR,ONE,V(istart,m),1,Lambda(1,m),1) ! Add (H0xK_i)^T * U_i CALL HessTR_Vec ( Hes0, U(istart,m), K(istart), Tmp ) CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1,m),1) END DO ! m=1:NADJ END DO ! Add H * dJac_dT_0^T * \sum(gamma_i U_i) ! Tmp holds sum gamma_i U_i !$$ IF (.NOT.Autonomous) THEN !$$ DO m = 1,NADJ !$$ Tmp(1:NVAR) = ZERO !$$ DO istage = 1, ros_S !$$ istart = NVAR*(istage-1) + 1 !$$ CALL WAXPY(NVAR,ros_Gamma(istage),U(istart,m),1,Tmp,1) !$$ END DO !$$#ifdef FULL_ALGEBRA !$$ Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp) !$$#else !$$ CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2) !$$#endif !$$ CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1,m),1) !$$ END DO ! m=1:NADJ !$$ END IF ! .NOT.Autonomous END DO TimeLoop !~~~> Save last state !~~~> Succesful exit IERR = 1 !~~~> The integration was successful !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ END SUBROUTINE ros_DadjRateInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_DadjInt ( & NADJ, Lambda, & Tstart, Tend, T, & !~~~> Error indicator IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the implementation of a generic RosenbrockADJ method ! defined by ros_S (no of stages) ! and its coefficients ros_{A,C,M,E,Alpha,Gamma} !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input: the initial condition at Tstart; Output: the solution at T INTEGER, INTENT(IN) :: NADJ !~~~> First order adjoint REAL(kind=dp), INTENT(INOUT) :: Lambda(NVAR,NADJ) !!~~~> Input: integration interval REAL(kind=dp), INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) REAL(kind=dp), INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables REAL(kind=dp) :: Ystage(NVAR*ros_S), K(NVAR*ros_S) REAL(kind=dp) :: U(NVAR*ros_S,NADJ), V(NVAR*ros_S,NADJ) #ifdef FULL_ALGEBRA REAL(kind=dp), DIMENSION(NVAR,NVAR) :: Jac, dJdT, Ghimj #else REAL(kind=dp), DIMENSION(LU_NONZERO) :: Jac, dJdT, Ghimj #endif REAL(kind=dp) :: Hes0(NHESS) REAL(kind=dp) :: Tmp(NVAR), Tmp2(NVAR) REAL(kind=dp) :: H, HC, HA, Tau INTEGER :: Pivot(NVAR), Direction INTEGER :: i, j, m, istage, istart, jstart !~~~> Local parameters REAL(kind=dp), PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 REAL(kind=dp), PARAMETER :: DeltaMin = 1.0d-5 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF (Tend >= Tstart) THEN Direction = +1 ELSE Direction = -1 END IF !~~~> Time loop begins below TimeLoop: DO WHILE ( stack_ptr > 0 ) !~~~> Recover checkpoints for stage values and vectors CALL ros_DPop( ros_S, T, H, Ystage, K, Ghimj, Pivot ) ! ISTATUS(Nstp) = ISTATUS(Nstp) + 1 !~~~> Compute LU decomposition IF (.NOT.SaveLU) THEN CALL JacTemplate(T,Ystage(1),Ghimj) ISTATUS(Njac) = ISTATUS(Njac) + 1 Tau = ONE/(Direction*H*ros_Gamma(1)) #ifdef FULL_ALGEBRA Ghimj(1:NVAR,1:NVAR) = -Ghimj(1:NVAR,1:NVAR) DO i=1,NVAR Ghimj(i,i) = Ghimj(i,i)+Tau END DO #else CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau END DO #endif CALL ros_Decomp( Ghimj, Pivot, j ) END IF !~~~> Compute Hessian at the beginning of the interval CALL HessTemplate(T,Ystage(1),Hes0) !~~~> Compute the stages Stage: DO istage = ros_S, 1, -1 !~~~> Current istage first entry istart = NVAR*(istage-1) + 1 !~~~> Compute U DO m = 1,NADJ CALL WCOPY(NVAR,Lambda(1,m),1,U(istart,m),1) CALL WSCAL(NVAR,ros_M(istage),U(istart,m),1) END DO ! m=1:NADJ DO j = istage+1, ros_S jstart = NVAR*(j-1) + 1 HA = ros_A((j-1)*(j-2)/2+istage) HC = ros_C((j-1)*(j-2)/2+istage)/(Direction*H) DO m = 1,NADJ CALL WAXPY(NVAR,HA,V(jstart,m),1,U(istart,m),1) CALL WAXPY(NVAR,HC,U(jstart,m),1,U(istart,m),1) END DO ! m=1:NADJ END DO DO m = 1,NADJ CALL ros_Solve('T', Ghimj, Pivot, U(istart,m)) END DO ! m=1:NADJ !~~~> Compute V Tau = T + ros_Alpha(istage)*Direction*H CALL JacTemplate(Tau,Ystage(istart),Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 DO m = 1,NADJ #ifdef FULL_ALGEBRA V(istart:istart+NVAR-1,m) = MATMUL(TRANSPOSE(Jac),U(istart:istart+NVAR-1,m)) #else CALL JacTR_SP_Vec(Jac,U(istart,m),V(istart,m)) #endif END DO ! m=1:NADJ END DO Stage IF (.NOT.Autonomous) THEN !~~~> Compute the Jacobian derivative with respect to T. ! Last "Jac" computed for stage 1 CALL ros_JacTimeDerivative ( T, Roundoff, Ystage(1), Jac, dJdT ) END IF !~~~> Compute the new solution !~~~> Compute Lambda DO istage=1,ros_S istart = NVAR*(istage-1) + 1 DO m = 1,NADJ ! Add V_i CALL WAXPY(NVAR,ONE,V(istart,m),1,Lambda(1,m),1) ! Add (H0xK_i)^T * U_i CALL HessTR_Vec ( Hes0, U(istart,m), K(istart), Tmp ) CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1,m),1) END DO ! m=1:NADJ END DO ! Add H * dJac_dT_0^T * \sum(gamma_i U_i) ! Tmp holds sum gamma_i U_i IF (.NOT.Autonomous) THEN DO m = 1,NADJ Tmp(1:NVAR) = ZERO DO istage = 1, ros_S istart = NVAR*(istage-1) + 1 CALL WAXPY(NVAR,ros_Gamma(istage),U(istart,m),1,Tmp,1) END DO #ifdef FULL_ALGEBRA Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp) #else CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2) #endif CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1,m),1) END DO ! m=1:NADJ END IF ! .NOT.Autonomous END DO TimeLoop !~~~> Save last state !~~~> Succesful exit IERR = 1 !~~~> The integration was successful !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ END SUBROUTINE ros_DadjInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_CadjInt ( & NADJ, Y, & Tstart, Tend, T, & AbsTol_adj, RelTol_adj, & !~~~> Error indicator IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the implementation of a generic RosenbrockADJ method ! defined by ros_S (no of stages) ! and its coefficients ros_{A,C,M,E,Alpha,Gamma} !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input: the initial condition at Tstart; Output: the solution at T INTEGER, INTENT(IN) :: NADJ REAL(kind=dp), INTENT(INOUT) :: Y(NVAR,NADJ) !~~~> Input: integration interval REAL(kind=dp), INTENT(IN) :: Tstart,Tend !~~~> Input: adjoint tolerances REAL(kind=dp), INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ) !~~~> Output: time at which the solution is returned (T=Tend if success) REAL(kind=dp), INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables REAL(kind=dp) :: Y0(NVAR) REAL(kind=dp) :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ) REAL(kind=dp) :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ) #ifdef FULL_ALGEBRA REAL(kind=dp), DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT #else REAL(kind=dp), DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT #endif REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau REAL(kind=dp) :: Err, Yerr(NVAR,NADJ) INTEGER :: Pivot(NVAR), Direction, ioffset, j, istage, iadj LOGICAL :: RejectLastH, RejectMoreH, Singular !~~~> Local parameters REAL(kind=dp), PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 REAL(kind=dp), PARAMETER :: DeltaMin = 1.0d-5 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Initial preparations T = Tstart RSTATUS(Nhexit) = 0.0_dp H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) ) IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin IF (Tend >= Tstart) THEN Direction = +1 ELSE Direction = -1 END IF H = Direction*H RejectLastH=.FALSE. RejectMoreH=.FALSE. !~~~> Time loop begins below TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff*abs(Tend) <= ZERO) & .OR. (Direction < 0).AND.((Tend-T)+Roundoff*abs(Tend) <= ZERO) ) ! Added *abs(Tend) by KS, A.Sandu for boundary cases IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps CALL ros_ErrorMsg(-6,T,H,IERR) RETURN END IF IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small CALL ros_ErrorMsg(-7,T,H,IERR) RETURN END IF !~~~> Limit H if necessary to avoid going beyond Tend RSTATUS(Nhexit) = H H = MIN(H,ABS(Tend-T)) !~~~> Interpolate forward solution CALL ros_cadj_Y( T, Y0 ) !~~~> Compute the Jacobian at current time CALL JacTemplate(T, Y0, Jac0) ISTATUS(Njac) = ISTATUS(Njac) + 1 !~~~> Compute the function derivative with respect to T IF (.NOT.Autonomous) THEN CALL ros_JacTimeDerivative ( T, Roundoff, Y0, & Jac0, dJdT ) DO iadj = 1, NADJ #ifdef FULL_ALGEBRA dFdT(1:NVAR,iadj) = MATMUL(TRANSPOSE(dJdT),Y(1:NVAR,iadj)) #else CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj)) #endif CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1) END DO END IF !~~~> Ydot = -J^T*Y #ifdef FULL_ALGEBRA Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) #else CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA Fcn0(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac0),Y(1:NVAR,iadj)) #else CALL JacTR_SP_Vec(Jac0,Y(1,iadj),Fcn0(1,iadj)) #endif END DO !~~~> Repeat step calculation until current step accepted UntilAccepted: DO CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), & Jac0,Ghimj,Pivot,Singular) IF (Singular) THEN ! More than 5 consecutive failed decompositions CALL ros_ErrorMsg(-8,T,H,IERR) RETURN END IF !~~~> Compute the stages Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR) ioffset = NVAR*(istage-1) ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN DO iadj = 1, NADJ CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1) END DO ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1) DO j = 1, istage-1 DO iadj = 1, NADJ CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) END DO END DO Tau = T + ros_Alpha(istage)*Direction*H ! CALL FunTemplate(Tau,Ynew,Fcn) ! ISTATUS(Nfun) = ISTATUS(Nfun) + 1 CALL ros_cadj_Y( Tau, Y0 ) CALL JacTemplate(Tau, Y0, Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR) #else CALL WSCAL(LU_NONZERO,(-ONE),Jac,1) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA Fcn(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac),Ynew(1:NVAR,iadj)) #else CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj)) #endif !CALL WSCAL(NVAR,(-ONE),Fcn(1,iadj),1) END DO END IF ! if istage == 1 elseif ros_NewF(istage) DO iadj = 1, NADJ CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1) END DO DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) DO iadj = 1, NADJ CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, & K(ioffset+1,iadj),1) END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) DO iadj = 1, NADJ CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1) END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj)) END DO END DO Stage !~~~> Compute the new solution DO iadj = 1, NADJ CALL WCOPY(NVAR,Y(1,iadj),1,Ynew(1,iadj),1) DO j=1,ros_S CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) END DO END DO !~~~> Compute the error estimation CALL WSCAL(NVAR*NADJ,ZERO,Yerr,1) DO j=1,ros_S DO iadj = 1, NADJ CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1,iadj),1,Yerr(1,iadj),1) END DO END DO !~~~> Max error among all adjoint components iadj = 1 Err = ros_ErrorNorm ( Y(1,iadj), Ynew(1,iadj), Yerr(1,iadj), & AbsTol_adj(1,iadj), RelTol_adj(1,iadj), VectorTol ) !~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO))) Hnew = H*Fac !~~~> Check the error magnitude and adjust step size ! ISTATUS(Nstp) = ISTATUS(Nstp) + 1 IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step ISTATUS(Nacc) = ISTATUS(Nacc) + 1 CALL WCOPY(NVAR*NADJ,Ynew,1,Y,1) T = T + Direction*H Hnew = MAX(Hmin,MIN(Hnew,Hmax)) IF (RejectLastH) THEN ! No step size increase after a rejected step Hnew = MIN(Hnew,H) END IF RSTATUS(Nhexit) = H RSTATUS(Nhnew) = Hnew RSTATUS(Ntexit) = T RejectLastH = .FALSE. RejectMoreH = .FALSE. H = Hnew EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED ELSE !~~~> Reject step IF (RejectMoreH) THEN Hnew = H*FacRej END IF RejectMoreH = RejectLastH RejectLastH = .TRUE. H = Hnew IF (ISTATUS(Nacc) >= 1) THEN ISTATUS(Nrej) = ISTATUS(Nrej) + 1 END IF END IF ! Err <= 1 END DO UntilAccepted END DO TimeLoop !~~~> Succesful exit IERR = 1 !~~~> The integration was successful END SUBROUTINE ros_CadjInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_SimpleCadjInt ( & NADJ, Y, & Tstart, Tend, T, & !~~~> Error indicator IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the implementation of a generic RosenbrockADJ method ! defined by ros_S (no of stages) ! and its coefficients ros_{A,C,M,E,Alpha,Gamma} !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input: the initial condition at Tstart; Output: the solution at T INTEGER, INTENT(IN) :: NADJ REAL(kind=dp), INTENT(INOUT) :: Y(NVAR,NADJ) !~~~> Input: integration interval REAL(kind=dp), INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) REAL(kind=dp), INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables REAL(kind=dp) :: Y0(NVAR) REAL(kind=dp) :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ) REAL(kind=dp) :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ) #ifdef FULL_ALGEBRA REAL(kind=dp),DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT #else REAL(kind=dp),DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT #endif REAL(kind=dp) :: H, HC, HG, Tau REAL(kind=dp) :: ghinv INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage, iadj INTEGER :: istack !~~~> Local parameters REAL(kind=dp), PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 REAL(kind=dp), PARAMETER :: DeltaMin = 1.0d-5 !~~~> Locally called functions ! REAL(kind=dp) WLAMCH ! EXTERNAL WLAMCH !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> INITIAL PREPARATIONS IF (Tend >= Tstart) THEN Direction = -1 ELSE Direction = +1 END IF !~~~> Time loop begins below TimeLoop: DO istack = stack_ptr,2,-1 T = chk_T(istack) H = chk_H(istack-1) !CALL WCOPY(NVAR,chk_Y(1,istack),1,Y0,1) Y0(1:NVAR) = chk_Y(1:NVAR,istack) !~~~> Compute the Jacobian at current time CALL JacTemplate(T, Y0, Jac0) ISTATUS(Njac) = ISTATUS(Njac) + 1 !~~~> Compute the function derivative with respect to T IF (.NOT.Autonomous) THEN CALL ros_JacTimeDerivative ( T, Roundoff, Y0, & Jac0, dJdT ) DO iadj = 1, NADJ #ifdef FULL_ALGEBRA dFdT(1:NVAR,iadj) = MATMUL(TRANSPOSE(dJdT),Y(1:NVAR,iadj)) #else CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj)) #endif CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1) END DO END IF !~~~> Ydot = -J^T*Y #ifdef FULL_ALGEBRA Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) #else CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA Fcn0(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac0),Y(1:NVAR,iadj)) #else CALL JacTR_SP_Vec(Jac0,Y(1,iadj),Fcn0(1,iadj)) #endif END DO !~~~> Construct Ghimj = 1/(H*ham) - Jac0 ghinv = ONE/(Direction*H*ros_Gamma(1)) #ifdef FULL_ALGEBRA Ghimj(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR) DO i=1,NVAR Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv END DO #endif !~~~> Compute LU decomposition CALL ros_Decomp( Ghimj, Pivot, j ) IF (j /= 0) THEN CALL ros_ErrorMsg(-8,T,H,IERR) PRINT*,' The matrix is singular !' STOP END IF !~~~> Compute the stages Stage: DO istage = 1, ros_S ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR) ioffset = NVAR*(istage-1) ! For the 1st istage the function has been computed previously IF ( istage == 1 ) THEN DO iadj = 1, NADJ CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1) END DO ! istage>1 and a new function evaluation is needed at the current istage ELSEIF ( ros_NewF(istage) ) THEN CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1) DO j = 1, istage-1 DO iadj = 1, NADJ CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), & K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1) END DO END DO Tau = T + ros_Alpha(istage)*Direction*H CALL ros_Hermite3( chk_T(istack-1), chk_T(istack), Tau, & chk_Y(1:NVAR,istack-1), chk_Y(1:NVAR,istack), & chk_dY(1:NVAR,istack-1), chk_dY(1:NVAR,istack), Y0 ) CALL JacTemplate(Tau, Y0, Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR) #else CALL WSCAL(LU_NONZERO,(-ONE),Jac,1) #endif DO iadj = 1, NADJ #ifdef FULL_ALGEBRA Fcn(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac),Ynew(1:NVAR,iadj)) #else CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj)) #endif END DO END IF ! if istage == 1 elseif ros_NewF(istage) DO iadj = 1, NADJ CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1) END DO DO j = 1, istage-1 HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H) DO iadj = 1, NADJ CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, & K(ioffset+1,iadj),1) END DO END DO IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN HG = Direction*H*ros_Gamma(istage) DO iadj = 1, NADJ CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1) END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj)) END DO END DO Stage !~~~> Compute the new solution DO iadj = 1, NADJ DO j=1,ros_S CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Y(1,iadj),1) END DO END DO END DO TimeLoop !~~~> Succesful exit IERR = 1 !~~~> The integration was successful END SUBROUTINE ros_SimpleCadjInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL(kind=dp) FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, & AbsTol, RelTol, VectorTol ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Computes the "scaled norm" of the error vector Yerr !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE ! Input arguments REAL(kind=dp), INTENT(IN) :: Y(NVAR), Ynew(NVAR), & Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR) LOGICAL, INTENT(IN) :: VectorTol ! Local variables REAL(kind=dp) :: Err, Scale, Ymax INTEGER :: i Err = ZERO DO i=1,NVAR Ymax = MAX(ABS(Y(i)),ABS(Ynew(i))) IF (VectorTol) THEN Scale = AbsTol(i)+RelTol(i)*Ymax ELSE Scale = AbsTol(1)+RelTol(1)*Ymax END IF Err = Err+(Yerr(i)/Scale)**2 END DO Err = SQRT(Err/NVAR) ros_ErrorNorm = MAX(Err,1.0d-10) END FUNCTION ros_ErrorNorm !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> The time partial derivative of the function by finite differences !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input arguments REAL(kind=dp), INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR) !~~~> Output arguments REAL(kind=dp), INTENT(OUT) :: dFdT(NVAR) !~~~> Local variables REAL(kind=dp) :: Delta REAL(kind=dp), PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6 Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL FunTemplate(T+Delta,Y,dFdT) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1) CALL WSCAL(NVAR,(ONE/Delta),dFdT,1) END SUBROUTINE ros_FunTimeDerivative !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, & Jac0, dJdT ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> The time partial derivative of the Jacobian by finite differences !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Arguments REAL(kind=dp), INTENT(IN) :: T, Roundoff, Y(NVAR) #ifdef FULL_ALGEBRA REAL(kind=dp), INTENT(IN) :: Jac0(NVAR,NVAR) REAL(kind=dp), INTENT(OUT) :: dJdT(NVAR,NVAR) #else REAL(kind=dp), INTENT(IN) :: Jac0(LU_NONZERO) REAL(kind=dp), INTENT(OUT) :: dJdT(LU_NONZERO) #endif !~~~> Local variables REAL(kind=dp) :: Delta Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL JacTemplate(T+Delta,Y,dJdT) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA CALL WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1) CALL WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1) #else CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1) CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1) #endif END SUBROUTINE ros_JacTimeDerivative !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, & Jac0, Ghimj, Pivot, Singular ) ! --- --- --- --- --- --- --- --- --- --- --- --- --- ! Prepares the LHS matrix for stage calculations ! 1. Construct Ghimj = 1/(H*gam) - Jac0 ! "(Gamma H) Inverse Minus Jacobian" ! 2. Repeat LU decomposition of Ghimj until successful. ! -half the step size if LU decomposition fails and retry ! -exit after 5 consecutive fails ! --- --- --- --- --- --- --- --- --- --- --- --- --- IMPLICIT NONE !~~~> Input arguments #ifdef FULL_ALGEBRA REAL(kind=dp), INTENT(IN) :: Jac0(NVAR,NVAR) #else REAL(kind=dp), INTENT(IN) :: Jac0(LU_NONZERO) #endif REAL(kind=dp), INTENT(IN) :: gam INTEGER, INTENT(IN) :: Direction !~~~> Output arguments #ifdef FULL_ALGEBRA REAL(kind=dp), INTENT(OUT) :: Ghimj(NVAR,NVAR) #else REAL(kind=dp), INTENT(OUT) :: Ghimj(LU_NONZERO) #endif LOGICAL, INTENT(OUT) :: Singular INTEGER, INTENT(OUT) :: Pivot(NVAR) !~~~> Inout arguments REAL(kind=dp), INTENT(INOUT) :: H ! step size is decreased when LU fails !~~~> Local variables INTEGER :: i, ising, Nconsecutive REAL(kind=dp) :: ghinv REAL(kind=dp), PARAMETER :: ONE = 1.0_dp, HALF = 0.5_dp Nconsecutive = 0 Singular = .TRUE. DO WHILE (Singular) !~~~> Construct Ghimj = 1/(H*gam) - Jac0 #ifdef FULL_ALGEBRA CALL WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1) CALL WSCAL(NVAR*NVAR,(-ONE),Ghimj,1) ghinv = ONE/(Direction*H*gam) DO i=1,NVAR Ghimj(i,i) = Ghimj(i,i)+ghinv END DO #else CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1) CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1) ghinv = ONE/(Direction*H*gam) DO i=1,NVAR Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv END DO #endif !~~~> Compute LU decomposition CALL ros_Decomp( Ghimj, Pivot, ising ) IF (ising == 0) THEN !~~~> If successful done Singular = .FALSE. ELSE ! ising .ne. 0 !~~~> If unsuccessful half the step size; if 5 consecutive fails then return ISTATUS(Nsng) = ISTATUS(Nsng) + 1 Nconsecutive = Nconsecutive+1 Singular = .TRUE. PRINT*,'Warning: LU Decomposition returned ising = ',ising IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions H = H*HALF ELSE ! More than 5 consecutive failed decompositions RETURN END IF ! Nconsecutive END IF ! ising END DO ! WHILE Singular END SUBROUTINE ros_PrepareMatrix !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_Decomp( A, Pivot, ising ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the LU decomposition !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Inout variables #ifdef FULL_ALGEBRA REAL(kind=dp), INTENT(INOUT) :: A(NVAR,NVAR) #else REAL(kind=dp), INTENT(INOUT) :: A(LU_NONZERO) #endif !~~~> Output variables INTEGER, INTENT(OUT) :: Pivot(NVAR), ising #ifdef FULL_ALGEBRA CALL DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising ) #else CALL KppDecomp ( A, ising ) Pivot(1) = 1 #endif ISTATUS(Ndec) = ISTATUS(Ndec) + 1 END SUBROUTINE ros_Decomp !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_Solve( How, A, Pivot, b ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the forward/backward substitution (using pre-computed LU decomposition) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input variables CHARACTER, INTENT(IN) :: How #ifdef FULL_ALGEBRA REAL(kind=dp), INTENT(IN) :: A(NVAR,NVAR) #else REAL(kind=dp), INTENT(IN) :: A(LU_NONZERO) #endif INTEGER, INTENT(IN) :: Pivot(NVAR) !~~~> InOut variables REAL(kind=dp), INTENT(INOUT) :: b(NVAR) SELECT CASE (How) CASE ('N') #ifdef FULL_ALGEBRA CALL DGETRS( 'N', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 ) #else CALL KppSolve( A, b ) #endif CASE ('T') #ifdef FULL_ALGEBRA CALL DGETRS( 'T', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 ) #else CALL KppSolveTR( A, b, b ) #endif CASE DEFAULT PRINT*,'Error: unknown argument in ros_Solve: How=',How STOP END SELECT ISTATUS(Nsol) = ISTATUS(Nsol) + 1 END SUBROUTINE ros_Solve !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_cadj_Y( T, Y ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Finds the solution Y at T by interpolating the stored forward trajectory !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input variables REAL(kind=dp), INTENT(IN) :: T !~~~> Output variables REAL(kind=dp), INTENT(OUT) :: Y(NVAR) !~~~> Local variables INTEGER :: i REAL(kind=dp), PARAMETER :: ONE = 1.0d0 ! chk_H, chk_T, chk_Y, chk_dY, chk_d2Y IF( (T < chk_T(1)).OR.(T> chk_T(stack_ptr)) ) THEN PRINT*,'Cannot locate solution at T = ',T PRINT*,'Stored trajectory is between Tstart = ',chk_T(1) PRINT*,' and Tend = ',chk_T(stack_ptr) STOP END IF DO i = 1, stack_ptr-1 IF( (T>= chk_T(i)).AND.(T<= chk_T(i+1)) ) EXIT END DO ! IF (.FALSE.) THEN ! ! CALL ros_Hermite5( chk_T(i), chk_T(i+1), T, & ! chk_Y(1,i), chk_Y(1,i+1), & ! chk_dY(1,i), chk_dY(1,i+1), & ! chk_d2Y(1,i), chk_d2Y(1,i+1), Y ) ! ! ELSE CALL ros_Hermite3( chk_T(i), chk_T(i+1), T, & chk_Y(1:NVAR,i), chk_Y(1:NVAR,i+1), & chk_dY(1:NVAR,i), chk_dY(1:NVAR,i+1), & Y ) ! ! END IF END SUBROUTINE ros_cadj_Y !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_Hermite3( a, b, T, Ya, Yb, Ja, Jb, Y ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for Hermite interpolation of order 5 on the interval [a,b] ! P = c(1) + c(2)*(x-a) + ... + c(4)*(x-a)^3 ! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb] !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input variables REAL(kind=dp), INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR) REAL(kind=dp), INTENT(IN) :: Ja(NVAR), Jb(NVAR) !~~~> Output variables REAL(kind=dp), INTENT(OUT) :: Y(NVAR) !~~~> Local variables REAL(kind=dp) :: Tau, amb(3), C(NVAR,4) REAL(kind=dp), PARAMETER :: ZERO = 0.0d0 INTEGER :: i, j amb(1) = 1.0d0/(a-b) DO i=2,3 amb(i) = amb(i-1)*amb(1) END DO ! c(1) = ya; CALL WCOPY(NVAR,Ya,1,C(1,1),1) ! c(2) = ja; CALL WCOPY(NVAR,Ja,1,C(1,2),1) ! c(3) = 2/(a-b)*ja + 1/(a-b)*jb - 3/(a - b)^2*ya + 3/(a - b)^2*yb ; CALL WCOPY(NVAR,Ya,1,C(1,3),1) CALL WSCAL(NVAR,-3.0*amb(2),C(1,3),1) CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1,3),1) CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1,3),1) CALL WAXPY(NVAR,amb(1),Jb,1,C(1,3),1) ! c(4) = 1/(a-b)^2*ja + 1/(a-b)^2*jb - 2/(a-b)^3*ya + 2/(a-b)^3*yb ; CALL WCOPY(NVAR,Ya,1,C(1,4),1) CALL WSCAL(NVAR,-2.0*amb(3),C(1,4),1) CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1,4),1) CALL WAXPY(NVAR,amb(2),Ja,1,C(1,4),1) CALL WAXPY(NVAR,amb(2),Jb,1,C(1,4),1) Tau = T - a CALL WCOPY(NVAR,C(1,4),1,Y,1) CALL WSCAL(NVAR,Tau**3,Y,1) DO j = 3,1,-1 CALL WAXPY(NVAR,TAU**(j-1),C(1,j),1,Y,1) END DO END SUBROUTINE ros_Hermite3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE ros_Hermite5( a, b, T, Ya, Yb, Ja, Jb, Ha, Hb, Y ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for Hermite interpolation of order 5 on the interval [a,b] ! P = c(1) + c(2)*(x-a) + ... + c(6)*(x-a)^5 ! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb], P"[a,b] = [Ha,Hb] !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Input variables REAL(kind=dp), INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR) REAL(kind=dp), INTENT(IN) :: Ja(NVAR), Jb(NVAR), Ha(NVAR), Hb(NVAR) !~~~> Output variables REAL(kind=dp), INTENT(OUT) :: Y(NVAR) !~~~> Local variables REAL(kind=dp) :: Tau, amb(5), C(NVAR,6) REAL(kind=dp), PARAMETER :: ZERO = 0.0d0, HALF = 0.5d0 INTEGER :: i, j amb(1) = 1.0d0/(a-b) DO i=2,5 amb(i) = amb(i-1)*amb(1) END DO ! c(1) = ya; CALL WCOPY(NVAR,Ya,1,C(1,1),1) ! c(2) = ja; CALL WCOPY(NVAR,Ja,1,C(1,2),1) ! c(3) = ha/2; CALL WCOPY(NVAR,Ha,1,C(1,3),1) CALL WSCAL(NVAR,HALF,C(1,3),1) ! c(4) = 10*amb(3)*ya - 10*amb(3)*yb - 6*amb(2)*ja - 4*amb(2)*jb + 1.5*amb(1)*ha - 0.5*amb(1)*hb ; CALL WCOPY(NVAR,Ya,1,C(1,4),1) CALL WSCAL(NVAR,10.0*amb(3),C(1,4),1) CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1,4),1) CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1,4),1) CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1,4),1) CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1,4),1) CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1,4),1) ! c(5) = 15*amb(4)*ya - 15*amb(4)*yb - 8.*amb(3)*ja - 7*amb(3)*jb + 1.5*amb(2)*ha - 1*amb(2)*hb ; CALL WCOPY(NVAR,Ya,1,C(1,5),1) CALL WSCAL(NVAR, 15.0*amb(4),C(1,5),1) CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1,5),1) CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1,5),1) CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1,5),1) CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1,5),1) CALL WAXPY(NVAR,-amb(2),Hb,1,C(1,5),1) ! c(6) = 6*amb(5)*ya - 6*amb(5)*yb - 3.*amb(4)*ja - 3.*amb(4)*jb + 0.5*amb(3)*ha -0.5*amb(3)*hb ; CALL WCOPY(NVAR,Ya,1,C(1,6),1) CALL WSCAL(NVAR, 6.0*amb(5),C(1,6),1) CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1,6),1) CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1,6),1) CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1,6),1) CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1,6),1) CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1,6),1) Tau = T - a CALL WCOPY(NVAR,C(1,6),1,Y,1) DO j = 5,1,-1 CALL WSCAL(NVAR,Tau,Y,1) CALL WAXPY(NVAR,ONE,C(1,j),1,Y,1) END DO END SUBROUTINE ros_Hermite5 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Ros2 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! --- AN L-STABLE METHOD, 2 stages, order 2 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE DOUBLE PRECISION g g = 1.0d0 + 1.0d0/SQRT(2.0d0) rosMethod = RS2 !~~~> Name of the method ros_Name = 'ROS-2' !~~~> Number of stages ros_S = 2 !~~~> The coefficient matrices A and C are strictly lower triangular. ! The lower triangular (subdiagonal) elements are stored in row-wise order: ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. ! The general mapping formula is: ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) ros_A(1) = (1.d0)/g ros_C(1) = (-2.d0)/g !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) ros_NewF(1) = .TRUE. ros_NewF(2) = .TRUE. !~~~> M_i = Coefficients for new step solution ros_M(1)= (3.d0)/(2.d0*g) ros_M(2)= (1.d0)/(2.d0*g) ! E_i = Coefficients for error estimator ros_E(1) = 1.d0/(2.d0*g) ros_E(2) = 1.d0/(2.d0*g) !~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus one ros_ELO = 2.0d0 !~~~> Y_stage_i ~ Y( T + H*Alpha_i ) ros_Alpha(1) = 0.0d0 ros_Alpha(2) = 1.0d0 !~~~> Gamma_i = \sum_j gamma_{i,j} ros_Gamma(1) = g ros_Gamma(2) =-g END SUBROUTINE Ros2 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Ros3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE rosMethod = RS3 !~~~> Name of the method ros_Name = 'ROS-3' !~~~> Number of stages ros_S = 3 !~~~> The coefficient matrices A and C are strictly lower triangular. ! The lower triangular (subdiagonal) elements are stored in row-wise order: ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. ! The general mapping formula is: ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) ros_A(1)= 1.d0 ros_A(2)= 1.d0 ros_A(3)= 0.d0 ros_C(1) = -0.10156171083877702091975600115545d+01 ros_C(2) = 0.40759956452537699824805835358067d+01 ros_C(3) = 0.92076794298330791242156818474003d+01 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) ros_NewF(1) = .TRUE. ros_NewF(2) = .TRUE. ros_NewF(3) = .FALSE. !~~~> M_i = Coefficients for new step solution ros_M(1) = 0.1d+01 ros_M(2) = 0.61697947043828245592553615689730d+01 ros_M(3) = -0.42772256543218573326238373806514d+00 ! E_i = Coefficients for error estimator ros_E(1) = 0.5d+00 ros_E(2) = -0.29079558716805469821718236208017d+01 ros_E(3) = 0.22354069897811569627360909276199d+00 !~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 ros_ELO = 3.0d0 !~~~> Y_stage_i ~ Y( T + H*Alpha_i ) ros_Alpha(1)= 0.0d+00 ros_Alpha(2)= 0.43586652150845899941601945119356d+00 ros_Alpha(3)= 0.43586652150845899941601945119356d+00 !~~~> Gamma_i = \sum_j gamma_{i,j} ros_Gamma(1)= 0.43586652150845899941601945119356d+00 ros_Gamma(2)= 0.24291996454816804366592249683314d+00 ros_Gamma(3)= 0.21851380027664058511513169485832d+01 END SUBROUTINE Ros3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Ros4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES ! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3 ! ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, ! SPRINGER-VERLAG (1990) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE rosMethod = RS4 !~~~> Name of the method ros_Name = 'ROS-4' !~~~> Number of stages ros_S = 4 !~~~> The coefficient matrices A and C are strictly lower triangular. ! The lower triangular (subdiagonal) elements are stored in row-wise order: ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. ! The general mapping formula is: ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) ros_A(1) = 0.2000000000000000d+01 ros_A(2) = 0.1867943637803922d+01 ros_A(3) = 0.2344449711399156d+00 ros_A(4) = ros_A(2) ros_A(5) = ros_A(3) ros_A(6) = 0.0D0 ros_C(1) =-0.7137615036412310d+01 ros_C(2) = 0.2580708087951457d+01 ros_C(3) = 0.6515950076447975d+00 ros_C(4) =-0.2137148994382534d+01 ros_C(5) =-0.3214669691237626d+00 ros_C(6) =-0.6949742501781779d+00 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) ros_NewF(1) = .TRUE. ros_NewF(2) = .TRUE. ros_NewF(3) = .TRUE. ros_NewF(4) = .FALSE. !~~~> M_i = Coefficients for new step solution ros_M(1) = 0.2255570073418735d+01 ros_M(2) = 0.2870493262186792d+00 ros_M(3) = 0.4353179431840180d+00 ros_M(4) = 0.1093502252409163d+01 !~~~> E_i = Coefficients for error estimator ros_E(1) =-0.2815431932141155d+00 ros_E(2) =-0.7276199124938920d-01 ros_E(3) =-0.1082196201495311d+00 ros_E(4) =-0.1093502252409163d+01 !~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 ros_ELO = 4.0d0 !~~~> Y_stage_i ~ Y( T + H*Alpha_i ) ros_Alpha(1) = 0.D0 ros_Alpha(2) = 0.1145640000000000d+01 ros_Alpha(3) = 0.6552168638155900d+00 ros_Alpha(4) = ros_Alpha(3) !~~~> Gamma_i = \sum_j gamma_{i,j} ros_Gamma(1) = 0.5728200000000000d+00 ros_Gamma(2) =-0.1769193891319233d+01 ros_Gamma(3) = 0.7592633437920482d+00 ros_Gamma(4) =-0.1049021087100450d+00 END SUBROUTINE Ros4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! --- A STIFFLY-STABLE METHOD, 4 stages, order 3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE rosMethod = RD3 !~~~> Name of the method ros_Name = 'RODAS-3' !~~~> Number of stages ros_S = 4 !~~~> The coefficient matrices A and C are strictly lower triangular. ! The lower triangular (subdiagonal) elements are stored in row-wise order: ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. ! The general mapping formula is: ! A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) ros_A(1) = 0.0d+00 ros_A(2) = 2.0d+00 ros_A(3) = 0.0d+00 ros_A(4) = 2.0d+00 ros_A(5) = 0.0d+00 ros_A(6) = 1.0d+00 ros_C(1) = 4.0d+00 ros_C(2) = 1.0d+00 ros_C(3) =-1.0d+00 ros_C(4) = 1.0d+00 ros_C(5) =-1.0d+00 ros_C(6) =-(8.0d+00/3.0d+00) !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) ros_NewF(1) = .TRUE. ros_NewF(2) = .FALSE. ros_NewF(3) = .TRUE. ros_NewF(4) = .TRUE. !~~~> M_i = Coefficients for new step solution ros_M(1) = 2.0d+00 ros_M(2) = 0.0d+00 ros_M(3) = 1.0d+00 ros_M(4) = 1.0d+00 !~~~> E_i = Coefficients for error estimator ros_E(1) = 0.0d+00 ros_E(2) = 0.0d+00 ros_E(3) = 0.0d+00 ros_E(4) = 1.0d+00 !~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 ros_ELO = 3.0d+00 !~~~> Y_stage_i ~ Y( T + H*Alpha_i ) ros_Alpha(1) = 0.0d+00 ros_Alpha(2) = 0.0d+00 ros_Alpha(3) = 1.0d+00 ros_Alpha(4) = 1.0d+00 !~~~> Gamma_i = \sum_j gamma_{i,j} ros_Gamma(1) = 0.5d+00 ros_Gamma(2) = 1.5d+00 ros_Gamma(3) = 0.0d+00 ros_Gamma(4) = 0.0d+00 END SUBROUTINE Rodas3 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Rodas4 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES ! ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, ! SPRINGER-VERLAG (1996) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE rosMethod = RD4 !~~~> Name of the method ros_Name = 'RODAS-4' !~~~> Number of stages ros_S = 6 !~~~> Y_stage_i ~ Y( T + H*Alpha_i ) ros_Alpha(1) = 0.000d0 ros_Alpha(2) = 0.386d0 ros_Alpha(3) = 0.210d0 ros_Alpha(4) = 0.630d0 ros_Alpha(5) = 1.000d0 ros_Alpha(6) = 1.000d0 !~~~> Gamma_i = \sum_j gamma_{i,j} ros_Gamma(1) = 0.2500000000000000d+00 ros_Gamma(2) =-0.1043000000000000d+00 ros_Gamma(3) = 0.1035000000000000d+00 ros_Gamma(4) =-0.3620000000000023d-01 ros_Gamma(5) = 0.0d0 ros_Gamma(6) = 0.0d0 !~~~> The coefficient matrices A and C are strictly lower triangular. ! The lower triangular (subdiagonal) elements are stored in row-wise order: ! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc. ! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j ) ! C(i,j) = ros_C( (i-1)*(i-2)/2 + j ) ros_A(1) = 0.1544000000000000d+01 ros_A(2) = 0.9466785280815826d+00 ros_A(3) = 0.2557011698983284d+00 ros_A(4) = 0.3314825187068521d+01 ros_A(5) = 0.2896124015972201d+01 ros_A(6) = 0.9986419139977817d+00 ros_A(7) = 0.1221224509226641d+01 ros_A(8) = 0.6019134481288629d+01 ros_A(9) = 0.1253708332932087d+02 ros_A(10) =-0.6878860361058950d+00 ros_A(11) = ros_A(7) ros_A(12) = ros_A(8) ros_A(13) = ros_A(9) ros_A(14) = ros_A(10) ros_A(15) = 1.0d+00 ros_C(1) =-0.5668800000000000d+01 ros_C(2) =-0.2430093356833875d+01 ros_C(3) =-0.2063599157091915d+00 ros_C(4) =-0.1073529058151375d+00 ros_C(5) =-0.9594562251023355d+01 ros_C(6) =-0.2047028614809616d+02 ros_C(7) = 0.7496443313967647d+01 ros_C(8) =-0.1024680431464352d+02 ros_C(9) =-0.3399990352819905d+02 ros_C(10) = 0.1170890893206160d+02 ros_C(11) = 0.8083246795921522d+01 ros_C(12) =-0.7981132988064893d+01 ros_C(13) =-0.3152159432874371d+02 ros_C(14) = 0.1631930543123136d+02 ros_C(15) =-0.6058818238834054d+01 !~~~> M_i = Coefficients for new step solution ros_M(1) = ros_A(7) ros_M(2) = ros_A(8) ros_M(3) = ros_A(9) ros_M(4) = ros_A(10) ros_M(5) = 1.0d+00 ros_M(6) = 1.0d+00 !~~~> E_i = Coefficients for error estimator ros_E(1) = 0.0d+00 ros_E(2) = 0.0d+00 ros_E(3) = 0.0d+00 ros_E(4) = 0.0d+00 ros_E(5) = 0.0d+00 ros_E(6) = 1.0d+00 !~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE) ! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE) ros_NewF(1) = .TRUE. ros_NewF(2) = .TRUE. ros_NewF(3) = .TRUE. ros_NewF(4) = .TRUE. ros_NewF(5) = .TRUE. ros_NewF(6) = .TRUE. !~~~> ros_ELO = estimator of local order - the minimum between the ! main and the embedded scheme orders plus 1 ros_ELO = 4.0d0 END SUBROUTINE Rodas4 END SUBROUTINE RosenbrockADJ ! and its internal procedures !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE FunTemplate( T, Y, Ydot ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the ODE function call. ! Updates the rate coefficients (and possibly the fixed species) at each call !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Input variables REAL(kind=dp), INTENT(IN) :: T, Y(NVAR) !~~~> Output variables REAL(kind=dp), INTENT(OUT) :: Ydot(NVAR) !~~~> Local variables REAL(kind=dp) :: Told !!$ Told = TIME !!$ TIME = T !!$ CALL Update_SUN() !!$ CALL Update_RCONST() CALL Fun( Y, FIX, RCONST, Ydot ) !!$ TIME = Told END SUBROUTINE FunTemplate !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE JacTemplate( T, Y, Jcb ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the ODE Jacobian call. ! Updates the rate coefficients (and possibly the fixed species) at each call !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Input variables REAL(kind=dp) :: T, Y(NVAR) !~~~> Output variables #ifdef FULL_ALGEBRA REAL(kind=dp) :: JV(LU_NONZERO), Jcb(NVAR,NVAR) #else REAL(kind=dp) :: Jcb(LU_NONZERO) #endif !~~~> Local variables REAL(kind=dp) :: Told #ifdef FULL_ALGEBRA INTEGER :: i, j #endif !!$ Told = TIME !!$ TIME = T !!$ CALL Update_SUN() !!$ CALL Update_RCONST() #ifdef FULL_ALGEBRA CALL Jac_SP(Y, FIX, RCONST, JV) DO j=1,NVAR DO i=1,NVAR Jcb(i,j) = 0.0_dp END DO END DO DO i=1,LU_NONZERO Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i) END DO #else CALL Jac_SP( Y, FIX, RCONST, Jcb ) #endif !!$ TIME = Told END SUBROUTINE JacTemplate !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE HessTemplate( T, Y, Hes ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the ODE Hessian call. ! Updates the rate coefficients (and possibly the fixed species) at each call !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Input variables REAL(kind=dp), INTENT(IN) :: T, Y(NVAR) !~~~> Output variables REAL(kind=dp), INTENT(OUT) :: Hes(NHESS) !~~~> Local variables REAL(kind=dp) :: Told !!$ Told = TIME !!$ TIME = T !!$ CALL Update_SUN() !!$ CALL Update_RCONST() CALL Hessian( Y, FIX, RCONST, Hes ) !!$ TIME = Told END SUBROUTINE HessTemplate END MODULE gckpp_adj_Integrator ! End of INTEGRATE function ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~