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

preparation of data struc. 

for adaptive time stepping
and cleaning
parent d411c604
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/adaptive_timestep_mod.F90 0 → 100644
+ 130
− 0
View file @ b3629a86
!! 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
src/advance_field_mod.F90
+ 1
− 1
View file @ b3629a86
......@@ -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
......
src/basic_mod.F90
+ 7
− 2
View file @ b3629a86
......@@ -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
......
src/time_integration_mod.F90
+ 228
− 147
View file @ b3629a86
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
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