!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! ! 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 ! ! This implementation is part of KPP - the Kinetic PreProcessor ! !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! MODULE KPP_ROOT_Integrator USE KPP_ROOT_Precision USE KPP_ROOT_Parameters USE KPP_ROOT_Global USE KPP_ROOT_LinearAlgebra USE KPP_ROOT_Rates USE KPP_ROOT_Function USE KPP_ROOT_Jacobian USE KPP_ROOT_Hessian USE KPP_ROOT_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 CONTAINS ! Routines in the module KPP_ROOT_Integrator !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, & ATOL_adj, RTOL_adj, & ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE !~~~> Y - Concentrations KPP_REAL :: 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 KPP_REAL :: Lambda(NVAR,NADJ) KPP_REAL, INTENT(IN) :: TIN ! TIN - Start Time KPP_REAL, INTENT(IN) :: TOUT ! TOUT - End Time !~~~> Tolerances for adjoint calculations ! (used only for full continuous adjoint) KPP_REAL, INTENT(IN) :: ATOL_adj(NVAR,NADJ), RTOL_adj(NVAR,NADJ) !~~~> Optional input parameters and statistics 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) KPP_REAL :: 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, & 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,')' 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, & 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) ) ! = 1 : no adjoint ! = 2 : discrete adjoint ( of method ICNTRL(3) ) ! = 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 KPP_REAL, INTENT(INOUT) :: Y(NVAR) INTEGER, INTENT(IN) :: NADJ KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ) KPP_REAL, INTENT(IN) :: Tstart,Tend KPP_REAL, INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR) KPP_REAL, INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ) INTEGER, INTENT(IN) :: ICNTRL(20) KPP_REAL, INTENT(IN) :: RCNTRL(20) INTEGER, INTENT(INOUT) :: ISTATUS(20) KPP_REAL, 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 KPP_REAL :: 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 INTEGER, PARAMETER :: bufsize = 200000 INTEGER :: stack_ptr = 0 ! last written entry KPP_REAL, DIMENSION(:), POINTER :: chk_H, chk_T KPP_REAL, DIMENSION(:,:), POINTER :: chk_Y, chk_K, chk_J KPP_REAL, DIMENSION(:,:), POINTER :: chk_dY, chk_d2Y !~~~> Local variables KPP_REAL :: Roundoff, FacMin, FacMax, FacRej, FacSafe KPP_REAL :: Hmin, Hmax, Hstart KPP_REAL :: Texit INTEGER :: i, UplimTol, Max_no_steps INTEGER :: AdjointType, CadjMethod LOGICAL :: Autonomous, VectorTol, SaveLU !~~~> Parameters KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, 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 !~~~> If Forward integration failed return IF (IERR<0) RETURN !~~~> 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_DadjInt ( & NADJ, Lambda, & 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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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 KPP_REAL :: T, H, Ystage(NVAR*S), K(NVAR*S) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: E(NVAR,NVAR) #else KPP_REAL :: 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 KPP_REAL :: T, H, Ystage(NVAR*S), K(NVAR*S) INTEGER :: P(NVAR) #ifdef FULL_ALGEBRA KPP_REAL :: E(NVAR,NVAR) #else KPP_REAL :: 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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ KPP_REAL :: 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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ KPP_REAL :: 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 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ KPP_REAL, 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 KPP_REAL, INTENT(INOUT) :: Y(NVAR) !~~~> Input: integration interval KPP_REAL, INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) KPP_REAL, INTENT(OUT) :: T !~~~> Input: tolerances KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR) !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables KPP_REAL :: Ynew(NVAR), Fcn0(NVAR), Fcn(NVAR) KPP_REAL :: K(NVAR*ros_S), dFdT(NVAR) KPP_REAL, DIMENSION(:), POINTER :: Ystage #ifdef FULL_ALGEBRA KPP_REAL :: Jac0(NVAR,NVAR), Ghimj(NVAR,NVAR) #else KPP_REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO) #endif KPP_REAL :: H, Hnew, HC, HG, Fac, Tau KPP_REAL :: 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 <= ZERO) & .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) ) 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 Fun_Template(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 Jac_Template(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:NVAR*j),1,Ynew,1) END DO Tau = T + ros_Alpha(istage)*Direction*H CALL Fun_Template(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:ioffset+NVAR),1) Ystage(ioffset+1:ioffset+NVAR) = Ynew(1:NVAR) END IF !CALL WCOPY(NVAR,Fcn,1,K(ioffset+1:ioffset+NVAR),1) K(ioffset+1:ioffset+NVAR) = Fcn(1:NVAR) 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:NVAR*j),1,K(ioffset+1:ioffset+NVAR),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:ioffset+NVAR),1) END IF CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR)) 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:NVAR*j),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:NVAR*j),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:NVAR) ) #endif IF (.NOT. Autonomous) THEN CALL WAXPY(NVAR,ONE,dFdT,1,K(1:NVAR),1) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1:NVAR) ) 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 Fun_Template(T,Y,Fcn0) ISTATUS(Nfun) = ISTATUS(Nfun) + 1 CALL Jac_Template(T,Y,Jac0) ISTATUS(Njac) = ISTATUS(Njac) + 1 #ifdef FULL_ALGEBRA K = MATMUL(Jac0,Fcn0) #else CALL Jac_SP_Vec( Jac0, Fcn0, K(1:NVAR) ) #endif IF (.NOT. Autonomous) THEN CALL ros_FunTimeDerivative ( T, Roundoff, Y, & Fcn0, dFdT ) CALL WAXPY(NVAR,ONE,dFdT,1,K(1:NVAR),1) END IF CALL ros_CPush( T, H, Y, Fcn0, K(1:NVAR) ) !~~~> 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_DadjInt ( & NADJ, Lambda, & Tstart, Tend, T, & !~~~> Error indicator IERR ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the implementation of a generic RosenbrockSOA 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 KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ) !!~~~> Input: integration interval KPP_REAL, INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) KPP_REAL, INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables KPP_REAL :: Ystage(NVAR*ros_S), K(NVAR*ros_S) KPP_REAL :: U(NVAR*ros_S,NADJ), V(NVAR*ros_S,NADJ) #ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: Jac, dJdT, Ghimj #else KPP_REAL, DIMENSION(LU_NONZERO) :: Jac, dJdT, Ghimj #endif KPP_REAL :: Hes0(NHESS) KPP_REAL :: Tmp(NVAR), Tmp2(NVAR) KPP_REAL :: H, HC, HA, Tau INTEGER :: Pivot(NVAR), Direction INTEGER :: i, j, m, istage, ioffset, joffset !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, 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 Jac_Template(T,Ystage(1:NVAR),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 Hess_Template(T,Ystage(1:NVAR),Hes0) !~~~> Compute the stages Stage: DO istage = ros_S, 1, -1 !~~~> Current istage first entry ioffset = NVAR*(istage-1) !~~~> Compute U DO m = 1,NADJ !CALL WCOPY(NVAR,Lambda(1:NVAR,m),1,U(ioffset+1:ioffset+NVAR,m),1) U(ioffset+1:ioffset+NVAR,m) = Lambda(1:NVAR,m) CALL WSCAL(NVAR,ros_M(istage),U(ioffset+1:ioffset+NVAR,m),1) END DO ! m=1:NADJ DO j = istage+1, ros_S joffset = NVAR*(j-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(joffset+1:joffset+NVAR,m),1,U(ioffset+1:ioffset+NVAR,m),1) CALL WAXPY(NVAR,HC,U(joffset+1:joffset+NVAR,m),1,U(ioffset+1:ioffset+NVAR,m),1) END DO ! m=1:NADJ END DO DO m = 1,NADJ CALL ros_Solve('T', Ghimj, Pivot, U(ioffset+1:ioffset+NVAR,m)) END DO ! m=1:NADJ !~~~> Compute V Tau = T + ros_Alpha(istage)*Direction*H CALL Jac_Template(Tau,Ystage(ioffset+1:ioffset+NVAR),Jac) ISTATUS(Njac) = ISTATUS(Njac) + 1 DO m = 1,NADJ #ifdef FULL_ALGEBRA V(ioffset+1:ioffset+NVAR,m) = MATMUL(TRANSPOSE(Jac),U(ioffset+1:ioffset+NVAR,m)) #else CALL JacTR_SP_Vec(Jac,U(ioffset+1:ioffset+NVAR,m),V(ioffset+1:ioffset+NVAR,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 ioffset = NVAR*(istage-1) DO m = 1,NADJ ! Add V_i CALL WAXPY(NVAR,ONE,V(ioffset+1:ioffset+NVAR,m),1,Lambda(1,m),1) ! Add (H0xK_i)^T * U_i CALL HessTR_Vec ( Hes0, U(ioffset+1:ioffset+NVAR,m), & K(ioffset+1:ioffset+NVAR), Tmp ) CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1:NVAR,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 ioffset = NVAR*(istage-1) CALL WAXPY(NVAR,ros_Gamma(istage),U(ioffset+1:ioffset+NVAR,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:NVAR,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 KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ) !~~~> Input: integration interval KPP_REAL, INTENT(IN) :: Tstart,Tend !~~~> Input: adjoint tolerances KPP_REAL, INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ) !~~~> Output: time at which the solution is returned (T=Tend if success) KPP_REAL, INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables KPP_REAL :: Y0(NVAR) KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ) KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ) #ifdef FULL_ALGEBRA KPP_REAL, DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT #else KPP_REAL, DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT #endif KPP_REAL :: H, Hnew, HC, HG, Fac, Tau KPP_REAL :: Err, Yerr(NVAR,NADJ) INTEGER :: Pivot(NVAR), Direction, ioffset, j, istage, iadj LOGICAL :: RejectLastH, RejectMoreH, Singular !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, 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 <= ZERO) & .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) ) 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 Jac_Template(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:NVAR,iadj),dFdT(1:NVAR,iadj)) #endif CALL WSCAL(NVAR,(-ONE),dFdT(1:NVAR,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:NVAR,iadj),Fcn0(1:NVAR,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) Fcn(1:NVAR,iadj) = Fcn0(1:NVAR,iadj) 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:NVAR*j,iadj),1,Ynew(1,iadj),1) END DO END DO Tau = T + ros_Alpha(istage)*Direction*H ! CALL Fun_Template(Tau,Ynew,Fcn) ! ISTATUS(Nfun) = ISTATUS(Nfun) + 1 CALL ros_cadj_Y( Tau, Y0 ) CALL Jac_Template(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:NVAR,iadj),Fcn(1:NVAR,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) K(ioffset+1:ioffset+NVAR,iadj) = Fcn(1:NVAR,iadj) 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:NVAR*j,iadj),1, & K(ioffset+1:ioffset+NVAR,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:NVAR,iadj),1,K(ioffset+1:ioffset+NVAR,iadj),1) END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR,iadj)) END DO END DO Stage !~~~> Compute the new solution DO iadj = 1, NADJ !CALL WCOPY(NVAR,Y(1:NVAR,iadj),1,Ynew(1,iadj),1) Ynew(1:NVAR,iadj) = Y(1:NVAR,iadj) DO j=1,ros_S CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1:NVAR*j,iadj),1,Ynew(1:NVAR,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:NVAR*j,iadj),1,Yerr(1:NVAR,iadj),1) END DO END DO !~~~> Max error among all adjoint components iadj = 1 Err = ros_ErrorNorm ( Y(1:NVAR,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 KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ) !~~~> Input: integration interval KPP_REAL, INTENT(IN) :: Tstart,Tend !~~~> Output: time at which the solution is returned (T=Tend if success) KPP_REAL, INTENT(OUT) :: T !~~~> Output: Error indicator INTEGER, INTENT(OUT) :: IERR ! ~~~~ Local variables KPP_REAL :: Y0(NVAR) KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ) KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ) #ifdef FULL_ALGEBRA KPP_REAL,DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT #else KPP_REAL,DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT #endif KPP_REAL :: H, HC, HG, Tau KPP_REAL :: ghinv INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage, iadj INTEGER :: istack !~~~> Local parameters KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0 KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5 !~~~> Locally called functions ! KPP_REAL 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 Jac_Template(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:NVAR,iadj),dFdT(1:NVAR,iadj)) #endif CALL WSCAL(NVAR,(-ONE),dFdT(1:NVAR,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:NVAR,iadj),Fcn0(1:NVAR,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:NVAR*j,iadj),1,Ynew(1:NVAR,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 Jac_Template(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:NVAR,iadj),Fcn(1:NVAR,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:NVAR*j,iadj),1, & K(ioffset+1:ioffset+NVAR,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:NVAR,iadj),1,K(ioffset+1:ioffset+NVAR,iadj),1) END DO END IF DO iadj = 1, NADJ CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR,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:NVAR*j,iadj),1,Y(1:NVAR,iadj),1) END DO END DO END DO TimeLoop !~~~> Succesful exit IERR = 1 !~~~> The integration was successful END SUBROUTINE ros_SimpleCadjInt !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, & AbsTol, RelTol, VectorTol ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Computes the "scaled norm" of the error vector Yerr !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPLICIT NONE ! Input arguments KPP_REAL, INTENT(IN) :: Y(NVAR), Ynew(NVAR), & Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR) LOGICAL, INTENT(IN) :: VectorTol ! Local variables KPP_REAL :: 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 KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR) !~~~> Output arguments KPP_REAL, INTENT(OUT) :: dFdT(NVAR) !~~~> Local variables KPP_REAL :: Delta KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6 Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL Fun_Template(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 KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR) #ifdef FULL_ALGEBRA KPP_REAL, INTENT(IN) :: Jac0(NVAR,NVAR) KPP_REAL, INTENT(OUT) :: dJdT(NVAR,NVAR) #else KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO) KPP_REAL, INTENT(OUT) :: dJdT(LU_NONZERO) #endif !~~~> Local variables KPP_REAL :: Delta Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T)) CALL Jac_Template(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 KPP_REAL, INTENT(IN) :: Jac0(NVAR,NVAR) #else KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO) #endif KPP_REAL, INTENT(IN) :: gam INTEGER, INTENT(IN) :: Direction !~~~> Output arguments #ifdef FULL_ALGEBRA KPP_REAL, INTENT(OUT) :: Ghimj(NVAR,NVAR) #else KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO) #endif LOGICAL, INTENT(OUT) :: Singular INTEGER, INTENT(OUT) :: Pivot(NVAR) !~~~> Inout arguments KPP_REAL, INTENT(INOUT) :: H ! step size is decreased when LU fails !~~~> Local variables INTEGER :: i, ising, Nconsecutive KPP_REAL :: ghinv KPP_REAL, 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 KPP_REAL, INTENT(INOUT) :: A(NVAR,NVAR) #else KPP_REAL, 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 KPP_REAL, INTENT(IN) :: A(NVAR,NVAR) #else KPP_REAL, INTENT(IN) :: A(LU_NONZERO) #endif INTEGER, INTENT(IN) :: Pivot(NVAR) !~~~> InOut variables KPP_REAL, 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 KPP_REAL, INTENT(IN) :: T !~~~> Output variables KPP_REAL, INTENT(OUT) :: Y(NVAR) !~~~> Local variables INTEGER :: i KPP_REAL, 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 KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR) KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR) !~~~> Output variables KPP_REAL, INTENT(OUT) :: Y(NVAR) !~~~> Local variables KPP_REAL :: Tau, amb(3), C(NVAR,4) KPP_REAL, 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:NVAR,1),1) ! c(2) = ja; CALL WCOPY(NVAR,Ja,1,C(1:NVAR,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:NVAR,3),1) CALL WSCAL(NVAR,-3.0*amb(2),C(1:NVAR,3),1) CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1:NVAR,3),1) CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1:NVAR,3),1) CALL WAXPY(NVAR,amb(1),Jb,1,C(1:NVAR,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:NVAR,4),1) CALL WSCAL(NVAR,-2.0*amb(3),C(1:NVAR,4),1) CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1:NVAR,4),1) CALL WAXPY(NVAR,amb(2),Ja,1,C(1:NVAR,4),1) CALL WAXPY(NVAR,amb(2),Jb,1,C(1:NVAR,4),1) Tau = T - a CALL WCOPY(NVAR,C(1:NVAR,4),1,Y,1) CALL WSCAL(NVAR,Tau**3,Y,1) DO j = 3,1,-1 CALL WAXPY(NVAR,TAU**(j-1),C(1:NVAR,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 KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR) KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR), Ha(NVAR), Hb(NVAR) !~~~> Output variables KPP_REAL, INTENT(OUT) :: Y(NVAR) !~~~> Local variables KPP_REAL :: Tau, amb(5), C(NVAR,6) KPP_REAL, 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:NVAR,1),1) ! c(2) = ja; CALL WCOPY(NVAR,Ja,1,C(1:NVAR,2),1) ! c(3) = ha/2; CALL WCOPY(NVAR,Ha,1,C(1:NVAR,3),1) CALL WSCAL(NVAR,HALF,C(1:NVAR,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:NVAR,4),1) CALL WSCAL(NVAR,10.0*amb(3),C(1:NVAR,4),1) CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1:NVAR,4),1) CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1:NVAR,4),1) CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1:NVAR,4),1) CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1:NVAR,4),1) CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1:NVAR,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:NVAR,5),1) CALL WSCAL(NVAR, 15.0*amb(4),C(1:NVAR,5),1) CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1:NVAR,5),1) CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1:NVAR,5),1) CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1:NVAR,5),1) CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1:NVAR,5),1) CALL WAXPY(NVAR,-amb(2),Hb,1,C(1:NVAR,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:NVAR,6),1) CALL WSCAL(NVAR, 6.0*amb(5),C(1:NVAR,6),1) CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1:NVAR,6),1) CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1:NVAR,6),1) CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1:NVAR,6),1) CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1:NVAR,6),1) CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1:NVAR,6),1) Tau = T - a CALL WCOPY(NVAR,C(1:NVAR,6),1,Y,1) DO j = 5,1,-1 CALL WSCAL(NVAR,Tau,Y,1) CALL WAXPY(NVAR,ONE,C(1:NVAR,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 Fun_Template( T, Y, Ydot ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the ODE function call. ! Updates the rate coefficients (and possibly the fixed species) at each call !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Input variables KPP_REAL, INTENT(IN) :: T, Y(NVAR) !~~~> Output variables KPP_REAL, INTENT(OUT) :: Ydot(NVAR) !~~~> Local variables KPP_REAL :: Told Told = TIME TIME = T CALL Update_SUN() CALL Update_RCONST() CALL Fun( Y, FIX, RCONST, Ydot ) TIME = Told END SUBROUTINE Fun_Template !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Jac_Template( T, Y, Jcb ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the ODE Jacobian call. ! Updates the rate coefficients (and possibly the fixed species) at each call !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Input variables KPP_REAL :: T, Y(NVAR) !~~~> Output variables #ifdef FULL_ALGEBRA KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR) #else KPP_REAL :: Jcb(LU_NONZERO) #endif !~~~> Local variables KPP_REAL :: 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 Jac_Template !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ SUBROUTINE Hess_Template( T, Y, Hes ) !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Template for the ODE Hessian call. ! Updates the rate coefficients (and possibly the fixed species) at each call !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !~~~> Input variables KPP_REAL, INTENT(IN) :: T, Y(NVAR) !~~~> Output variables KPP_REAL, INTENT(OUT) :: Hes(NHESS) !~~~> Local variables KPP_REAL :: Told Told = TIME TIME = T CALL Update_SUN() CALL Update_RCONST() CALL Hessian( Y, FIX, RCONST, Hes ) TIME = Told END SUBROUTINE Hess_Template END MODULE KPP_ROOT_Integrator