Skip to content
Snippets Groups Projects
Commit 13ea6942 authored by Antoine Cyril David Hoffmann's avatar Antoine Cyril David Hoffmann :seedling:
Browse files

Added multiple time stepping schemes

parent 33cad120
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/time_integration_mod.F90
+ 153
− 119
View file @ 13ea6942
...@@ -60,56 +60,104 @@ CONTAINS ...@@ -60,56 +60,104 @@ CONTAINS
IMPLICIT NONE IMPLICIT NONE
SELECT CASE (numerical_scheme) SELECT CASE (numerical_scheme)
CASE ('RK4') ! Order 2 methods
CALL RK4 CASE ('RK2')
CASE ('SSPx_RK3') CALL RK2
CALL SSPx_RK3
CASE ('SSPx_RK2') CASE ('SSPx_RK2')
CALL SSPx_RK2 CALL SSPx_RK2
! Order 3 methods
CASE ('RK3')
CALL RK3
CASE ('SSP_RK3')
CALL SSP_RK3
CASE ('SSPx_RK3')
CALL SSPx_RK3
CASE ('IMEX_SSP2')
CALL IMEX_SSP2
CASE ('ARK2')
CALL ARK2
! Order 4 methods
CASE ('RK4')
CALL RK4
! Order 5 methods
CASE ('DOPRI5') CASE ('DOPRI5')
CALL DOPRI5 CALL DOPRI5
CASE DEFAULT CASE DEFAULT
IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.' IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
END SELECT END SELECT
IF (my_id .EQ. 0) WRITE(*,*) " Time integration with ", numerical_scheme IF (my_id .EQ. 0) WRITE(*,*) " Time integration with ", numerical_scheme
END SUBROUTINE set_numerical_scheme END SUBROUTINE set_numerical_scheme
SUBROUTINE RK4 !!! second order time schemes
! Butcher coeff for RK4 (default) SUBROUTINE RK2
! Butcher coeff for clasical RK2 (Heun's)
USE basic USE basic
USE prec_const USE prec_const
IMPLICIT NONE IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 4 INTEGER,PARAMETER :: nbstep = 2
CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 2
ntimelevel = 4
c_E(1) = 0.0_dp c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp c_E(2) = 1.0_dp
c_E(3) = 1.0_dp/2.0_dp b_E(1) = 1._dp/2._dp
c_E(4) = 1.0_dp b_E(2) = 1._dp/2._dp
A_E(2,1) = 1._dp
END SUBROUTINE RK2
b_E(1) = 1.0_dp/6.0_dp SUBROUTINE SSPx_RK2
b_E(2) = 1.0_dp/3.0_dp ! DOESNT WORK
b_E(3) = 1.0_dp/3.0_dp ! Butcher coeff for modified strong stability preserving RK2
b_E(4) = 1.0_dp/6.0_dp ! used in GX (Hammett 2022, Mandell et al. 2022)
USE basic
USE prec_const
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 2
REAL(dp) :: alpha, beta
alpha = 1._dp/SQRT(2._dp)
beta = SQRT(2._dp) - 1._dp
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 2
c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp
b_E(1) = alpha*beta/2._dp
b_E(2) = alpha/2._dp
A_E(2,1) = alpha
! b_E(1) = 1._dp
! b_E(2) = 1._dp/SQRT(2._dp)
! A_E(2,1) = 1._dp/SQRT(2._dp)
END SUBROUTINE SSPx_RK2
!!! third order time schemes
SUBROUTINE RK3
! Butcher coeff for classical RK3
USE basic
USE prec_const
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp
c_E(3) = 1.0_dp
b_E(1) = 1._dp/6._dp
b_E(2) = 2._dp/3._dp
b_E(3) = 1._dp/6._dp
A_E(2,1) = 1.0_dp/2.0_dp A_E(2,1) = 1.0_dp/2.0_dp
A_E(3,2) = 1.0_dp/2.0_dp A_E(3,1) = -1._dp
A_E(4,3) = 1.0_dp A_E(3,2) = 2._dp
END SUBROUTINE RK3
END SUBROUTINE RK4
SUBROUTINE SSPx_RK3 SUBROUTINE SSPx_RK3
! DOESNT WORK
! Butcher coeff for modified strong stability preserving RK3 ! Butcher coeff for modified strong stability preserving RK3
! used in GX (Hammett 2022, Mandell et al. 2022) ! used in GX (Hammett 2022, Mandell et al. 2022)
USE basic USE basic
USE prec_const USE prec_const
IMPLICIT NONE IMPLICIT NONE
...@@ -125,125 +173,122 @@ CONTAINS ...@@ -125,125 +173,122 @@ CONTAINS
! w2 = 0.3769892220587804931852815570891834795475_dp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2 ! w2 = 0.3769892220587804931852815570891834795475_dp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
w3 = 1._dp/a1 - w2 * (1._dp + w1) w3 = 1._dp/a1 - w2 * (1._dp + w1)
! w3 = 1.3368459739528868457369981115334667265415_dp ! w3 = 1.3368459739528868457369981115334667265415_dp
CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 3 ntimelevel = 3
c_E(1) = 0.0_dp c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp c_E(2) = 1.0_dp/2.0_dp
c_E(3) = 1.0_dp/2.0_dp c_E(3) = 1.0_dp/2.0_dp
b_E(1) = a1 * (w1*w2 + w3) b_E(1) = a1 * (w1*w2 + w3)
b_E(2) = a2 * w2 b_E(2) = a2 * w2
b_E(3) = a3 b_E(3) = a3
A_E(2,1) = a1 A_E(2,1) = a1
A_E(3,1) = a1 * w1 A_E(3,1) = a1 * w1
A_E(3,2) = a2 A_E(3,2) = a2
END SUBROUTINE SSPx_RK3 END SUBROUTINE SSPx_RK3
SUBROUTINE SSPx_RK2 SUBROUTINE IMEX_SSP2
! Butcher coeff for modified strong stability preserving RK2 !! Version of Rokhzadi 2017 (An Optimally Stable and Accurate Second-Order
! used in GX (Hammett 2022, Mandell et al. 2022) ! SSP Runge-Kutta IMEX Scheme for Atmospheric Applications)
USE basic USE basic
USE prec_const USE prec_const
IMPLICIT NONE IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 2 INTEGER,PARAMETER :: nbstep = 3
REAL(dp) :: alpha, beta
alpha = 1._dp/SQRT(2._dp)
beta = SQRT(2._dp) - 1._dp
CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
ntimelevel = 2 c_E(1) = 0._dp
c_E(2) = 0.711664700366941_dp
c_E(1) = 0.0_dp c_E(3) = 0.994611536833690_dp
c_E(2) = 1.0_dp/2.0_dp b_E(1) = 0.398930808264688_dp
b_E(2) = 0.345755244189623_dp
b_E(1) = alpha*beta/2._dp b_E(3) = 0.255313947545689_dp
b_E(2) = alpha/2._dp A_E(2,1) = 0.711664700366941_dp
A_E(3,1) = 0.077338168947683_dp
A_E(2,1) = alpha A_E(3,2) = 0.917273367886007_dp
END SUBROUTINE IMEX_SSP2
END SUBROUTINE SSPx_RK2
SUBROUTINE ARK2
SUBROUTINE DOPRI5 !! Version of Rokhzadi 2017 (An Optimally Stable and Accurate Second-Order
! Butcher coeff for DOPRI5 --> Stiffness detection ! SSP Runge-Kutta IMEX Scheme for Atmospheric Applications)
! DOPRI5 used for stiffness detection.
! 5 order method/7 stages
USE basic USE basic
USE prec_const
IMPLICIT NONE IMPLICIT NONE
INTEGER,PARAMETER :: nbstep =7 INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
c_E(1) = 0._dp
c_E(2) = 2._dp*(1._dp - 1._dp/SQRT2)
c_E(3) = 1._dp
b_E(1) = 1._dp/(2._dp*SQRT2)
b_E(2) = 1._dp/(2._dp*SQRT2)
b_E(3) = 1._dp - 1._dp/SQRT2
A_E(2,1) = 2._dp*(1._dp - 1._dp/SQRT2)
A_E(3,1) = 1._dp - (3._dp + 2._dp*SQRT2)/6._dp
A_E(3,2) = (3._dp + 2._dp*SQRT2)/6._dp
END SUBROUTINE ARK2
SUBROUTINE SSP_RK3
! Butcher coeff for strong stability preserving RK3
USE basic
USE prec_const
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
c_E(1) = 0.0_dp
c_E(2) = 1.0_dp
c_E(3) = 1.0_dp/2.0_dp
b_E(1) = 1._dp/6._dp
b_E(2) = 1._dp/6._dp
b_E(3) = 2._dp/3._dp
A_E(2,1) = 1._dp
A_E(3,1) = 1._dp/4._dp
A_E(3,2) = 1._dp/4._dp
END SUBROUTINE SSP_RK3
!!! fourth order time schemes
SUBROUTINE RK4
! Butcher coeff for RK4 (default)
USE basic
USE prec_const
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 4
CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 4
c_E(1) = 0.0_dp
c_E(2) = 1.0_dp/2.0_dp
c_E(3) = 1.0_dp/2.0_dp
c_E(4) = 1.0_dp
b_E(1) = 1.0_dp/6.0_dp
b_E(2) = 1.0_dp/3.0_dp
b_E(3) = 1.0_dp/3.0_dp
b_E(4) = 1.0_dp/6.0_dp
A_E(2,1) = 1.0_dp/2.0_dp
A_E(3,2) = 1.0_dp/2.0_dp
A_E(4,3) = 1.0_dp
END SUBROUTINE RK4
ntimelevel = 7 !!! fifth order time schemes
! SUBROUTINE DOPRI5
c_E(1) = 0._dp
c_E(2) = 1.0_dp/5.0_dp
c_E(3) = 3.0_dp /10.0_dp
c_E(4) = 4.0_dp/5.0_dp
c_E(5) = 8.0_dp/9.0_dp
c_E(6) = 1.0_dp
c_E(7) = 1.0_dp
!
A_E(2,1) = 1.0_dp/5.0_dp
A_E(3,1) = 3.0_dp/40.0_dp
A_E(3,2) = 9.0_dp/40.0_dp
A_E(4,1) = 44.0_dp/45.0_dp
A_E(4,2) = -56.0_dp/15.0_dp
A_E(4,3) = 32.0_dp/9.0_dp
A_E(5,1 ) = 19372.0_dp/6561.0_dp
A_E(5,2) = -25360.0_dp/2187.0_dp
A_E(5,3) = 64448.0_dp/6561.0_dp
A_E(5,4) = -212.0_dp/729.0_dp
A_E(6,1) = 9017.0_dp/3168.0_dp
A_E(6,2)= -355.0_dp/33.0_dp
A_E(6,3) = 46732.0_dp/5247.0_dp
A_E(6,4) = 49.0_dp/176.0_dp
A_E(6,5) = -5103.0_dp/18656.0_dp
A_E(7,1) = 35.0_dp/384.0_dp
A_E(7,3) = 500.0_dp/1113.0_dp
A_E(7,4) = 125.0_dp/192.0_dp
A_E(7,5) = -2187.0_dp/6784.0_dp
A_E(7,6) = 11.0_dp/84.0_dp
!
b_E(1) = 35.0_dp/384.0_dp
b_E(2) = 0._dp
b_E(3) = 500.0_dp/1113.0_dp
b_E(4) = 125.0_dp/192.0_dp
b_E(5) = -2187.0_dp/6784.0_dp
b_E(6) = 11.0_dp/84.0_dp
b_E(7) = 0._dp
!
END SUBROUTINE DOPRI5
SUBROUTINE DOPRI5_ADAPT
! Butcher coeff for DOPRI5 --> Stiffness detection ! Butcher coeff for DOPRI5 --> Stiffness detection
! DOPRI5 used for stiffness detection. ! DOPRI5 used for stiffness detection.
! 5 order method/7 stages ! 5 order method/7 stages
USE basic USE basic
IMPLICIT NONE IMPLICIT NONE
INTEGER,PARAMETER :: nbstep =7 INTEGER,PARAMETER :: nbstep =7
CALL allocate_array_dp1(c_E,1,nbstep) CALL allocate_array_dp1(c_E,1,nbstep)
CALL allocate_array_dp1(b_E,1,nbstep) CALL allocate_array_dp1(b_E,1,nbstep)
CALL allocate_array_dp1(b_Es,1,nbstep)
CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep) CALL allocate_array_dp2(A_E,1,nbstep,1,nbstep)
ntimelevel = 7 ntimelevel = 7
!
c_E(1) = 0._dp c_E(1) = 0._dp
c_E(2) = 1.0_dp/5.0_dp c_E(2) = 1.0_dp/5.0_dp
c_E(3) = 3.0_dp /10.0_dp c_E(3) = 3.0_dp /10.0_dp
...@@ -251,7 +296,6 @@ CONTAINS ...@@ -251,7 +296,6 @@ CONTAINS
c_E(5) = 8.0_dp/9.0_dp c_E(5) = 8.0_dp/9.0_dp
c_E(6) = 1.0_dp c_E(6) = 1.0_dp
c_E(7) = 1.0_dp c_E(7) = 1.0_dp
!
A_E(2,1) = 1.0_dp/5.0_dp A_E(2,1) = 1.0_dp/5.0_dp
A_E(3,1) = 3.0_dp/40.0_dp A_E(3,1) = 3.0_dp/40.0_dp
A_E(3,2) = 9.0_dp/40.0_dp A_E(3,2) = 9.0_dp/40.0_dp
...@@ -272,7 +316,6 @@ CONTAINS ...@@ -272,7 +316,6 @@ CONTAINS
A_E(7,4) = 125.0_dp/192.0_dp A_E(7,4) = 125.0_dp/192.0_dp
A_E(7,5) = -2187.0_dp/6784.0_dp A_E(7,5) = -2187.0_dp/6784.0_dp
A_E(7,6) = 11.0_dp/84.0_dp A_E(7,6) = 11.0_dp/84.0_dp
!
b_E(1) = 35.0_dp/384.0_dp b_E(1) = 35.0_dp/384.0_dp
b_E(2) = 0._dp b_E(2) = 0._dp
b_E(3) = 500.0_dp/1113.0_dp b_E(3) = 500.0_dp/1113.0_dp
...@@ -280,15 +323,6 @@ CONTAINS ...@@ -280,15 +323,6 @@ CONTAINS
b_E(5) = -2187.0_dp/6784.0_dp b_E(5) = -2187.0_dp/6784.0_dp
b_E(6) = 11.0_dp/84.0_dp b_E(6) = 11.0_dp/84.0_dp
b_E(7) = 0._dp b_E(7) = 0._dp
! END SUBROUTINE DOPRI5
b_Es(1) = 5179.0_dp/57600.0_dp
b_Es(2) = 0._dp
b_Es(3) = 7571.0_dp/16695.0_dp
b_Es(4) = 393.0_dp/640.0_dp
b_Es(5) = -92097.0_dp/339200.0_dp
b_Es(6) = 187.0_dp/2100.0_dp
b_Es(7) = 1._dp/40._dp
!
END SUBROUTINE DOPRI5_ADAPT
END MODULE time_integration END MODULE time_integration
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment