From b3629a86bd452ed7b046b9e33e54e50e179adaaf Mon Sep 17 00:00:00 2001
From: Antoine Hoffmann <antoine.hoffmann@epfl.ch>
Date: Tue, 4 Jul 2023 15:54:48 +0200
Subject: [PATCH] preparation of data struc. for adaptive time stepping and
cleaning
---
src/adaptive_timestep_mod.F90 | 130 ++++++++++++
src/advance_field_mod.F90 | 2 +-
src/basic_mod.F90 | 9 +-
src/time_integration_mod.F90 | 375 +++++++++++++++++++++-------------
4 files changed, 366 insertions(+), 150 deletions(-)
create mode 100644 src/adaptive_timestep_mod.F90
diff --git a/src/adaptive_timestep_mod.F90 b/src/adaptive_timestep_mod.F90
new file mode 100644
index 00000000..c8581ea6
--- /dev/null
+++ b/src/adaptive_timestep_mod.F90
@@ -0,0 +1,130 @@
+!! This module is a compilation of the routines implemented in GBS by L. Stenger
+!! for adaptive time stepping.
+MODULE adaptive_timestep
+ USE prec_const, ONLY: xp
+
+ IMPLICIT NONE
+
+ logical, public, protected :: time_scheme_is_adaptive = .false.
+ logical, public, protected :: should_discard_step = .false.
+ real(xp) :: adaptive_accuracy_ratio = 1._dp
+ integer :: adaptive_error_estimator_order = -1
+ real(xp) :: dt_min = 1e-8_xp !< Minimum allowable timestep for adaptive time schemes
+ real(xp) :: dt_max = 1e1_xp !< Maximum allowable timestep for adaptive time schemes
+ real(xp) :: adaptive_safety = 0.9_xp !< Safety factor of the adaptive timestep controller
+ !> @{
+ !> @brief Error tolerance for adaptive time schemes
+ !> @details
+ !> Typically, an adaptive time scheme will provide an error estimation for a
+ !> given time step "witdh". The allowable error will be defined as
+ !> `atol + rtol * abs(f)` (where `f` the quantity which is being integrated).
+ real(xp), public, protected :: adaptive_error_atol=1e-6_dp, adaptive_error_rtol=1e-3_dp
+ !> @}
+ integer, public, protected :: ntimelevel_rhs=4 !< number of time levels required for the rhs
+
+
+ CONTAINS
+
+
+ SUBROUTINE adaptive_timestep_readinputs
+ ! Read the input parameters
+ USE prec_const
+ USE basic, ONLY : lu_in
+ IMPLICIT NONE
+ namelist /TIME_INTEGRATION_PAR/ adaptive_safety, adaptive_error_atol, adaptive_error_rtol
+
+ READ(lu_in,time_integration_par)
+ CALL set_numerical_scheme
+ END SUBROUTINE adaptive_timestep_readinputs
+
+ !> Computes the next time step for adaptive time schemes
+ !>
+ !> Provided a relative stepping error, defined by `tolerance/error` where
+ !> `tolerance` is the admissible error and `error` the estimated error of the
+ !> current step, compute an optimal `dt` for the next step.
+ !>
+ !> @param[in] dt Current time step
+ !> @param[in] errscale Normalized error of the current step
+ !> @returns Optimal `dt` for the next step
+ pure function adaptive_compute_next_dt(dt, errscale) result(next_dt)
+ real(xp), intent(in) :: dt, errscale
+ real(xp), parameter :: MINSCALE=.1_xp, MAXSCALE=5._xp
+ real(xp) :: next_dt, error_exponent
+ ! The parameters in this routine seem to work "fine". Attempts to be too
+ ! ambitious, especially with `adaptive_safety`, usually increase the rate
+ ! at which time steps are discarded!
+ if(errscale > 1._dp) then
+ error_exponent = 1 / real(adaptive_error_estimator_order + 1, xp)
+ else
+ error_exponent = 1 / real(adaptive_error_estimator_order, xp)
+ end if
+ next_dt = dt * max(MINSCALE, min(MAXSCALE, adaptive_safety * errscale**error_exponent))
+ end function adaptive_compute_next_dt
+
+ !> Setup for adaptive time schemes
+ !>
+ !> Prepares the next full embedded Runge-Kutta step (not substep) by adjusting
+ !> basic::dt based on the previous time step, along with extra cleanup in case
+ !> the previous step was discarded.
+ subroutine adaptive_time_scheme_setup()
+ use basic, only: cstepP, dt, meW, stepP
+ real(xp) :: old_dt
+
+ old_dt = dt
+ call set_dt(adaptive_compute_next_dt(dt, adaptive_accuracy_ratio))
+
+ if(should_discard_step) then
+ ! Note this relates to the previous step!
+ if(meW==0 .and. dt/=old_dt) print *, "step discarded. dt update ", old_dt, " => ", dt
+ if(meW==0 .and. dt/=old_dt) print *, "errscale was ", adaptive_accuracy_ratio
+ end if
+
+ ! Setting `adaptive_accuracy_ratio` is equivalent to saying that the step
+ ! has zero error, i.e. "assume the step is perfect now, and update later
+ ! when actually computing the error after performing the step"
+ should_discard_step = .false.
+ adaptive_accuracy_ratio = huge(adaptive_accuracy_ratio)
+ end subroutine adaptive_time_scheme_setup
+
+ !> Updates the normalized adaptive time scheme error
+ !>
+ !> This routine is useful to combine the normalized error of the full system
+ !> of equation by recording the worst (highest) error, which will eventually
+ !> be passed to time_integration::adaptive_compute_next_dt.
+ !>
+ !> @param[in] errscale Normalized time step error
+ !>
+ !> @sa time_integration::adaptive_compute_next_dt
+ subroutine adaptive_set_error(errscale)
+ real(xp), intent(in) :: errscale
+ logical :: satisfies_tolerance
+ satisfies_tolerance = errscale > 1._xp
+ adaptive_accuracy_ratio = min(adaptive_accuracy_ratio, errscale)
+ should_discard_step = should_discard_step .or. (.not. satisfies_tolerance)
+ end subroutine adaptive_set_error
+
+ !> Returns `.true.` if the current step should be discarded.
+ pure function adaptive_must_recompute_step() result(res)
+ logical :: res
+ res = should_discard_step
+ end function adaptive_must_recompute_step
+
+ subroutine set_updatetlevel(new_updatetlevel)
+ integer, intent(in) :: new_updatetlevel
+ updatetlevel = new_updatetlevel
+ end subroutine set_updatetlevel
+
+ subroutine set_updatetlevel_rhs
+ select case (numerical_scheme)
+ case ('rkc2')
+ if (updatetlevel == 1)then
+ updatetlevel_rhs = 1
+ else
+ updatetlevel_rhs = 2
+ end if
+ case default
+ updatetlevel_rhs = updatetlevel
+ end select
+ end subroutine set_updatetlevel_rhs
+
+END MODULE adaptive_timestep
\ No newline at end of file
diff --git a/src/advance_field_mod.F90 b/src/advance_field_mod.F90
index 9ee8c4be..d67d080c 100644
--- a/src/advance_field_mod.F90
+++ b/src/advance_field_mod.F90
@@ -36,7 +36,7 @@ CONTAINS
DO ia =1,local_na
IF( evolve_mom(ipi,iji) )&
moments(ia,ipi,iji,iky,ikx,izi,1) = moments(ia,ipi,iji,iky,ikx,izi,1) &
- + dt*b_E(istage)*moments_rhs(ia,ip,ij,iky,ikx,iz,istage)
+ + dt*b_E(istage,1)*moments_rhs(ia,ip,ij,iky,ikx,iz,istage)
END DO
END DO
END DO
diff --git a/src/basic_mod.F90 b/src/basic_mod.F90
index 704aad2a..0c45dd45 100644
--- a/src/basic_mod.F90
+++ b/src/basic_mod.F90
@@ -45,7 +45,7 @@ MODULE basic
! Routines interfaces
PUBLIC :: allocate_array, basic_outputinputs,basic_data,&
speak, str, increase_step, increase_cstep, increase_time, display_h_min_s,&
- set_basic_cp, daytim, start_chrono, stop_chrono
+ set_basic_cp, daytim, start_chrono, stop_chrono, change_dt
! Interface for allocating arrays, these routines allocate and initialize directly to zero
INTERFACE allocate_array
MODULE PROCEDURE allocate_array_xp1,allocate_array_xp2,allocate_array_xp3, &
@@ -123,6 +123,11 @@ CONTAINS
IMPLICIT NONE
cstep = cstep + 1
END SUBROUTINE
+ SUBROUTINE change_dt(new_dt)
+ IMPLICIT NONE
+ REAL(xp), INTENT(IN) :: new_dt
+ dt = new_dt
+ END SUBROUTINE
SUBROUTINE increase_time
IMPLICIT NONE
time = time + dt
@@ -244,7 +249,7 @@ CONTAINS
! "Convert an integer to string."
integer, intent(in) :: k
character(len=10) :: str_
- write (str_, "(i2.2)") k
+ write (str_, "(I8)") k
str_ = adjustl(str_)
end function str_int
diff --git a/src/time_integration_mod.F90 b/src/time_integration_mod.F90
index e1a2aa51..5b8e5654 100644
--- a/src/time_integration_mod.F90
+++ b/src/time_integration_mod.F90
@@ -1,33 +1,52 @@
MODULE time_integration
USE prec_const
+
IMPLICIT NONE
PRIVATE
INTEGER, PUBLIC, PROTECTED :: ntimelevel=4 ! Total number of time levels required by the numerical scheme
INTEGER, PUBLIC, PROTECTED :: updatetlevel ! Current time level to be updated
- real(xp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E,A_I
- real(xp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: b_E,b_Es,b_I
- real(xp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E,c_I !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
+ real(xp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: A_E
+ real(xp),PUBLIC,PROTECTED,DIMENSION(:,:),ALLOCATABLE :: b_E
+ real(xp),PUBLIC,PROTECTED,DIMENSION(:),ALLOCATABLE :: c_E !Coeff for Expl/Implic time integration in case of time dependent RHS (i.e. dy/dt = f(y,t)) see Baptiste Frei CSE Rapport 06/17
character(len=10),PUBLIC,PROTECTED :: numerical_scheme='RK4'
+ !!---- start: This part is taken from GBS adaptive time stepping by L. Stenger
+ ! For more information check gbs.epfl.ch
+ logical, public, protected :: time_scheme_is_adaptive = .false.
+ logical, public, protected :: should_discard_step = .false.
+ real(xp) :: adaptive_accuracy_ratio = 1._xp
+ integer :: adaptive_error_estimator_order = -1
+ real(xp) :: dt_min = 1e-8_xp !< Minimum allowable timestep for adaptive time schemes
+ real(xp) :: dt_max = 1e1_xp !< Maximum allowable timestep for adaptive time schemes
+ real(xp) :: adaptive_safety = 0.9_xp !< Safety factor of the adaptive timestep controller
+ !> @{
+ !> @brief Error tolerance for adaptive time schemes
+ !> @details
+ !> Typically, an adaptive time scheme will provide an error estimation for a
+ !> given time step "witdh". The allowable error will be defined as
+ !> `atol + rtol * abs(f)` (where `f` the quantity which is being integrated).
+ real(xp), public, protected :: adaptive_error_atol=1e-6_xp, adaptive_error_rtol=1e-3_xp
+ !> @}
+ integer, public, protected :: ntimelevel_rhs = 4 !< number of time levels required for the rhs
+ integer, public, protected :: updatetlevel_rhs = 1 !< time level to be updated for rhs
+ !!---- end
+
PUBLIC :: set_updatetlevel, time_integration_readinputs, time_integration_outputinputs
CONTAINS
- SUBROUTINE set_updatetlevel(new_updatetlevel)
- INTEGER, INTENT(in) :: new_updatetlevel
- updatetlevel = new_updatetlevel
- END SUBROUTINE set_updatetlevel
-
SUBROUTINE time_integration_readinputs
! Read the input parameters
USE prec_const
USE basic, ONLY : lu_in
IMPLICIT NONE
NAMELIST /TIME_INTEGRATION_PAR/ numerical_scheme
+ namelist /TIME_INTEGRATION_PAR/ adaptive_safety, adaptive_error_atol, adaptive_error_rtol
+
READ(lu_in,time_integration_par)
CALL set_numerical_scheme
END SUBROUTINE time_integration_readinputs
@@ -52,24 +71,21 @@ CONTAINS
! Order 2 methods
CASE ('RK2')
CALL RK2
- CASE ('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
+ ! Adaptative scheme
+ CASE ('ODE23')
+ time_scheme_is_adaptive = .true.
+ CALL ode23
! Order 4 methods
CASE ('RK4')
CALL RK4
! Order 5 methods
CASE ('DOPRI5')
+ time_scheme_is_adaptive = .true.
CALL DOPRI5
CASE DEFAULT
IF (my_id .EQ. 0) WRITE(*,*) 'Cannot initialize time integration scheme. Name invalid.'
@@ -85,41 +101,16 @@ CONTAINS
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 2
CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
+ CALL allocate_array(b_E,1,nbstep,1,1)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 2
- c_E(1) = 0.0_xp
- c_E(2) = 1.0_xp
- b_E(1) = 1._xp/2._xp
- b_E(2) = 1._xp/2._xp
+ c_E(1) = 0.0_xp
+ c_E(2) = 1.0_xp
+ b_E(1,1) = 1._xp/2._xp
+ b_E(2,1) = 1._xp/2._xp
A_E(2,1) = 1._xp
END SUBROUTINE RK2
- SUBROUTINE SSPx_RK2
- ! DOESNT WORK
- ! Butcher coeff for modified strong stability preserving RK2
- ! used in GX (Hammett 2022, Mandell et al. 2022)
- USE basic
- USE prec_const
- IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep = 2
- REAL(xp) :: alpha, beta
- alpha = 1._xp/SQRT(2._xp)
- beta = SQRT(2._xp) - 1._xp
- CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
- CALL allocate_array(A_E,1,nbstep,1,nbstep)
- ntimelevel = 2
- c_E(1) = 0.0_xp
- c_E(2) = 1.0_xp/2.0_xp
- b_E(1) = alpha*beta/2._xp
- b_E(2) = alpha/2._xp
- A_E(2,1) = alpha
- ! b_E(1) = 1._xp
- ! b_E(2) = 1._xp/SQRT(2._xp)
- ! A_E(2,1) = 1._xp/SQRT(2._xp)
- END SUBROUTINE SSPx_RK2
-
!!! third order time schemes
SUBROUTINE RK3
! Butcher coeff for classical RK3
@@ -128,98 +119,20 @@ CONTAINS
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
+ CALL allocate_array(b_E,1,nbstep,1,1)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
c_E(1) = 0.0_xp
c_E(2) = 1.0_xp/2.0_xp
c_E(3) = 1.0_xp
- b_E(1) = 1._xp/6._xp
- b_E(2) = 2._xp/3._xp
- b_E(3) = 1._xp/6._xp
+ b_E(1,1) = 1._xp/6._xp
+ b_E(2,1) = 2._xp/3._xp
+ b_E(3,1) = 1._xp/6._xp
A_E(2,1) = 1.0_xp/2.0_xp
A_E(3,1) = -1._xp
A_E(3,2) = 2._xp
END SUBROUTINE RK3
- SUBROUTINE SSPx_RK3
- ! DOESNT WORK
- ! Butcher coeff for modified strong stability preserving RK3
- ! used in GX (Hammett 2022, Mandell et al. 2022)
- USE basic
- USE prec_const
- IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep = 3
- REAL(xp) :: a1, a2, a3, w1, w2, w3
- a1 = (1._xp/6._xp)**(1._xp/3._xp)! (1/6)^(1/3)
- ! a1 = 0.5503212081491044571635029569733887910843_xp ! (1/6)^(1/3)
- a2 = a1
- a3 = a1
- w1 = 0.5_xp*(-1._xp + SQRT( 9._xp - 2._xp * 6._xp**(2._xp/3._xp))) ! (-1 + sqrt(9-2*6^(2/3)))/2
- ! w1 = 0.2739744023885328783052273138309828937054_xp ! (sqrt(9-2*6^(2/3))-1)/2
- w2 = 0.5_xp*(-1._xp + 6._xp**(2._xp/3._xp) - SQRT(9._xp - 2._xp * 6._xp**(2._xp/3._xp))) ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
- ! w2 = 0.3769892220587804931852815570891834795475_xp ! (6^(2/3)-1-sqrt(9-2*6^(2/3)))/2
- w3 = 1._xp/a1 - w2 * (1._xp + w1)
- ! w3 = 1.3368459739528868457369981115334667265415_xp
- CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
- CALL allocate_array(A_E,1,nbstep,1,nbstep)
- ntimelevel = 3
- c_E(1) = 0.0_xp
- c_E(2) = 1.0_xp/2.0_xp
- c_E(3) = 1.0_xp/2.0_xp
- b_E(1) = a1 * (w1*w2 + w3)
- b_E(2) = a2 * w2
- b_E(3) = a3
- A_E(2,1) = a1
- A_E(3,1) = a1 * w1
- A_E(3,2) = a2
- END SUBROUTINE SSPx_RK3
-
- SUBROUTINE IMEX_SSP2
- !! Version of Rokhzadi 2017 (An Optimally Stable and Accurate Second-Order
- ! SSP Runge-Kutta IMEX Scheme for Atmospheric Applications)
- USE basic
- USE prec_const
- IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep = 3
- CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
- CALL allocate_array(A_E,1,nbstep,1,nbstep)
- ntimelevel = 3
- c_E(1) = 0._xp
- c_E(2) = 0.711664700366941_xp
- c_E(3) = 0.994611536833690_xp
- b_E(1) = 0.398930808264688_xp
- b_E(2) = 0.345755244189623_xp
- b_E(3) = 0.255313947545689_xp
- A_E(2,1) = 0.711664700366941_xp
- A_E(3,1) = 0.077338168947683_xp
- A_E(3,2) = 0.917273367886007_xp
- END SUBROUTINE IMEX_SSP2
-
- SUBROUTINE ARK2
- !! Version of Rokhzadi 2017 (An Optimally Stable and Accurate Second-Order
- ! SSP Runge-Kutta IMEX Scheme for Atmospheric Applications)
- USE basic
- USE prec_const
- IMPLICIT NONE
- INTEGER,PARAMETER :: nbstep = 3
- CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
- CALL allocate_array(A_E,1,nbstep,1,nbstep)
- ntimelevel = 3
- c_E(1) = 0._xp
- c_E(2) = 2._xp*(1._xp - 1._xp/SQRT2)
- c_E(3) = 1._xp
- b_E(1) = 1._xp/(2._xp*SQRT2)
- b_E(2) = 1._xp/(2._xp*SQRT2)
- b_E(3) = 1._xp - 1._xp/SQRT2
- A_E(2,1) = 2._xp*(1._xp - 1._xp/SQRT2)
- A_E(3,1) = 1._xp - (3._xp + 2._xp*SQRT2)/6._xp
- A_E(3,2) = (3._xp + 2._xp*SQRT2)/6._xp
- END SUBROUTINE ARK2
-
SUBROUTINE SSP_RK3
! Butcher coeff for strong stability preserving RK3
USE basic
@@ -227,15 +140,15 @@ CONTAINS
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 3
CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
+ CALL allocate_array(b_E,1,nbstep,1,1)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 3
c_E(1) = 0.0_xp
c_E(2) = 1.0_xp
c_E(3) = 1.0_xp/2.0_xp
- b_E(1) = 1._xp/6._xp
- b_E(2) = 1._xp/6._xp
- b_E(3) = 2._xp/3._xp
+ b_E(1,1) = 1._xp/6._xp
+ b_E(2,1) = 1._xp/6._xp
+ b_E(3,1) = 2._xp/3._xp
A_E(2,1) = 1._xp
A_E(3,1) = 1._xp/4._xp
A_E(3,2) = 1._xp/4._xp
@@ -249,17 +162,17 @@ CONTAINS
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep = 4
CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
+ CALL allocate_array(b_E,1,nbstep,1,1)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
ntimelevel = 4
c_E(1) = 0.0_xp
c_E(2) = 1.0_xp/2.0_xp
c_E(3) = 1.0_xp/2.0_xp
c_E(4) = 1.0_xp
- b_E(1) = 1.0_xp/6.0_xp
- b_E(2) = 1.0_xp/3.0_xp
- b_E(3) = 1.0_xp/3.0_xp
- b_E(4) = 1.0_xp/6.0_xp
+ b_E(1,1) = 1.0_xp/6.0_xp
+ b_E(2,1) = 1.0_xp/3.0_xp
+ b_E(3,1) = 1.0_xp/3.0_xp
+ b_E(4,1) = 1.0_xp/6.0_xp
A_E(2,1) = 1.0_xp/2.0_xp
A_E(3,2) = 1.0_xp/2.0_xp
A_E(4,3) = 1.0_xp
@@ -274,20 +187,26 @@ CONTAINS
IMPLICIT NONE
INTEGER,PARAMETER :: nbstep =7
CALL allocate_array(c_E,1,nbstep)
- CALL allocate_array(b_E,1,nbstep)
+ CALL allocate_array(b_E,1,nbstep,1,2)
CALL allocate_array(A_E,1,nbstep,1,nbstep)
+
+ adaptive_error_estimator_order = 4
ntimelevel = 7
- c_E(1) = 0._xp
+ ntimelevel_rhs = 7
+
+ c_E = 0._xp
c_E(2) = 1.0_xp/5.0_xp
c_E(3) = 3.0_xp /10.0_xp
c_E(4) = 4.0_xp/5.0_xp
c_E(5) = 8.0_xp/9.0_xp
c_E(6) = 1.0_xp
c_E(7) = 1.0_xp
+
+ A_E = 0._xp
A_E(2,1) = 1.0_xp/5.0_xp
A_E(3,1) = 3.0_xp/40.0_xp
A_E(3,2) = 9.0_xp/40.0_xp
- A_E(4,1) = 44.0_xp/45.0_xp
+ A_E(4,1) = 44.0_xp/45.0_xp
A_E(4,2) = -56.0_xp/15.0_xp
A_E(4,3) = 32.0_xp/9.0_xp
A_E(5,1 ) = 19372.0_xp/6561.0_xp
@@ -304,13 +223,175 @@ CONTAINS
A_E(7,4) = 125.0_xp/192.0_xp
A_E(7,5) = -2187.0_xp/6784.0_xp
A_E(7,6) = 11.0_xp/84.0_xp
- b_E(1) = 35.0_xp/384.0_xp
- b_E(2) = 0._xp
- b_E(3) = 500.0_xp/1113.0_xp
- b_E(4) = 125.0_xp/192.0_xp
- b_E(5) = -2187.0_xp/6784.0_xp
- b_E(6) = 11.0_xp/84.0_xp
- b_E(7) = 0._xp
+
+ b_E = 0._xp
+ b_E(1, 1) = 35._xp/384._xp
+ b_E(2, 1) = 0._xp
+ b_E(3, 1) = 500._xp/1113._xp
+ b_E(4, 1) = 125._xp/192._xp
+ b_E(5, 1) = -2187._xp/6784._xp
+ b_E(6, 1) = 11._xp/84._xp
+ b_E(7, 1) = 0._xp
+ b_E(1, 2) = 5179._xp/57600._xp
+ b_E(2, 2) = 0._xp
+ b_E(3, 2) = 7571._xp/16695._xp
+ b_E(4, 2) = 393._xp/640._xp
+ b_E(5, 2) = -92097._xp/339200._xp
+ b_E(6, 2) = 187._xp/2100._xp
+ b_E(7, 2) = 1._xp/40._xp
END SUBROUTINE DOPRI5
+ !!-------------------------------------------------------------------------
+ !!---- This part is taken from GBS adaptive time stepping by L. Stenger ---
+ !> Updates the time step
+ !>
+ !> If the new time step in not within the closed interval `[dt_min, dt_max]`,
+ !> the simulation will be terminated
+ !>
+ !> @param[in] new_dt Desired time step value
+ subroutine set_dt(new_dt)
+ use basic, only: dt, nlend, change_dt
+ use parallel, ONLY: comm_p, my_id
+ real(xp), intent(in) :: new_dt
+ integer :: ierr
+ logical :: mlend
+ mlend = .false.
+ if(dt<=dt_min .or. dt>=dt_max) then
+ mlend = .true.
+ if(my_id == 0) then
+ print "(A, 3(G12.5,A))", "The adaptive time integration scheme tried to set the time step to ", dt, &
+ ", which is outside of the allowed range of [", dt_min, ", ", dt_max, &
+ "]." // new_line("A") // "The simulation will be terminated."
+ end if
+ else
+ CALL change_dt(new_dt)
+ end if
+ call mpi_allreduce(mlend, nlend, 1, MPI_LOGICAL, MPI_LOR, comm_p, ierr)
+ end subroutine set_dt
+ !> Computes the next time step for adaptive time schemes
+ !>
+ !> Provided a relative stepping error, defined by `tolerance/error` where
+ !> `tolerance` is the admissible error and `error` the estimated error of the
+ !> current step, compute an optimal `dt` for the next step.
+ !>
+ !> @param[in] dt Current time step
+ !> @param[in] errscale Normalized error of the current step
+ !> @returns Optimal `dt` for the next step
+ pure function adaptive_compute_next_dt(dt, errscale) result(next_dt)
+ real(xp), intent(in) :: dt, errscale
+ real(xp), parameter :: MINSCALE=.1_xp, MAXSCALE=5._xp
+ real(xp) :: next_dt, error_exponent
+ ! The parameters in this routine seem to work "fine". Attempts to be too
+ ! ambitious, especially with `adaptive_safety`, usually increase the rate
+ ! at which time steps are discarded!
+ if(errscale > 1._xp) then
+ error_exponent = 1 / real(adaptive_error_estimator_order + 1, xp)
+ else
+ error_exponent = 1 / real(adaptive_error_estimator_order, xp)
+ end if
+ next_dt = dt * max(MINSCALE, min(MAXSCALE, adaptive_safety * errscale**error_exponent))
+ end function adaptive_compute_next_dt
+
+ !> Setup for adaptive time schemes
+ !>
+ !> Prepares the next full embedded Runge-Kutta step (not substep) by adjusting
+ !> basic::dt based on the previous time step, along with extra cleanup in case
+ !> the previous step was discarded.
+ subroutine adaptive_time_scheme_setup()
+ use basic, only: dt, str, speak
+ real(xp) :: old_dt
+
+ old_dt = dt
+ call set_dt(adaptive_compute_next_dt(dt, adaptive_accuracy_ratio))
+
+ if(should_discard_step) then
+ ! Note this relates to the previous step!
+ if(dt/=old_dt) CALL speak("step discarded. dt update "//str(old_dt)//" => "//str(dt))
+ if(dt/=old_dt) CALL speak("errscale was "//str(adaptive_accuracy_ratio))
+ end if
+
+ ! Setting `adaptive_accuracy_ratio` is equivalent to saying that the step
+ ! has zero error, i.e. "assume the step is perfect now, and update later
+ ! when actually computing the error after performing the step"
+ should_discard_step = .false.
+ adaptive_accuracy_ratio = huge(adaptive_accuracy_ratio)
+ end subroutine adaptive_time_scheme_setup
+
+ !> Updates the normalized adaptive time scheme error
+ !>
+ !> This routine is useful to combine the normalized error of the full system
+ !> of equation by recording the worst (highest) error, which will eventually
+ !> be passed to time_integration::adaptive_compute_next_dt.
+ !>
+ !> @param[in] errscale Normalized time step error
+ !>
+ !> @sa time_integration::adaptive_compute_next_dt
+ subroutine adaptive_set_error(errscale)
+ real(xp), intent(in) :: errscale
+ logical :: satisfies_tolerance
+ satisfies_tolerance = errscale > 1._xp
+ adaptive_accuracy_ratio = min(adaptive_accuracy_ratio, errscale)
+ should_discard_step = should_discard_step .or. (.not. satisfies_tolerance)
+ end subroutine adaptive_set_error
+
+ !> Returns `.true.` if the current step should be discarded.
+ pure function adaptive_must_recompute_step() result(res)
+ logical :: res
+ res = should_discard_step
+ end function adaptive_must_recompute_step
+
+ subroutine set_updatetlevel(new_updatetlevel)
+ integer, intent(in) :: new_updatetlevel
+ updatetlevel = new_updatetlevel
+ end subroutine set_updatetlevel
+
+ subroutine set_updatetlevel_rhs
+ select case (numerical_scheme)
+ case ('rkc2')
+ if (updatetlevel == 1)then
+ updatetlevel_rhs = 1
+ else
+ updatetlevel_rhs = 2
+ end if
+ case default
+ updatetlevel_rhs = updatetlevel
+ end select
+ end subroutine set_updatetlevel_rhs
+
+ !> Sets up Butcher coefficients for the Bogacki-Shampine method (orders 3 and 2)
+ subroutine ode23
+ use basic
+ integer,parameter :: nbstep = 4
+ CALL allocate_array(c_E,1,nbstep)
+ CALL allocate_array(b_E,1,nbstep, 1, 2)
+ CALL allocate_array(A_E,1,nbstep,1,nbstep)
+
+ adaptive_error_estimator_order = 2
+ ntimelevel = 4
+ ntimelevel_rhs = 4
+
+ c_E = 0._xp
+ c_E(2) = 1.0_xp/2.0_xp
+ c_E(3) = 3.0_xp /4.0_xp
+ c_E(4) = 1._xp
+
+ A_E = 0._xp
+ A_E(2,1) = 1.0_xp/2.0_xp
+ A_E(3,2) = 3.0_xp/4.0_xp
+ A_E(4,1) = 2.0_xp/9.0_xp
+ A_E(4,2) = 1.0_xp/3.0_xp
+ A_E(4,3) = 4.0_xp/9.0_xp
+
+ b_E = 0._xp
+ b_E(1, 1) = 2._xp/9._xp
+ b_E(2, 1) = 1._xp/3._xp
+ b_E(3, 1) = 4._xp/9._xp
+ b_E(4, 1) = 0._xp
+ b_E(1, 2) = 7._xp/24._xp
+ b_E(2, 2) = 1._xp/4._xp
+ b_E(3, 2) = 1._xp/3._xp
+ b_E(4, 2) = 1._xp/8._xp
+ end subroutine ode23
+ !!-------------------------------------------------------------------------
+ !!-------------------------------------------------------------------------
END MODULE time_integration
--
GitLab