Sola/GEOS-Chem-adjoint-v35-note
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
2370ddd7939201be4feaa9287b28c2842edecabb
GEOS-Chem-adjoint-v35-note/code/adjoint/gckpp_adj_Integrator.f90
Xuesong (Steve) f842cc9639 Add files via upload
2018-08-28 00:33:48 -04:00

2906 lines
96 KiB
Fortran
Raw Blame History

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! 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
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Reference in New Issue View Git Blame Copy Permalink