EVOLUTION: Implement estimation of the initial time step.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2021-10-19 13:11:41 -03:00
parent 6ca36c93e3
commit d501b19041

View File

@ -70,6 +70,7 @@ module evolution
logical , save :: error_control = .false.
character(len=255), save :: name_int = ""
character(len=255), save :: name_norm = ""
integer , save :: order = 2
integer , save :: stages = 2
integer , save :: registers = 1
real(kind=8) , save :: cfl = 5.0d-01
@ -118,7 +119,7 @@ module evolution
! declare public subroutines
!
public :: initialize_evolution, finalize_evolution, print_evolution
public :: advance, new_time_step
public :: initialize_time_step, advance, new_time_step
! declare public variables
!
@ -222,12 +223,14 @@ module evolution
case ("euler", "EULER", "rk1", "RK1")
evolve => evolve_euler
order = 1
registers = 1
name_int = "1st order Euler"
case ("rk2", "RK2")
evolve => evolve_rk2
order = 2
registers = 2
name_int = "2nd order Runge-Kutta"
@ -235,6 +238,7 @@ module evolution
error_control = .true.
evolve => evolve_ssprk2_m
order = 2
registers = 3
stages = max(2, min(9, stages))
cfl = (stages - 1) * cfl
@ -243,12 +247,14 @@ module evolution
case ("rk3", "RK3")
evolve => evolve_rk3
order = 3
registers = 2
name_int = "3rd order Runge-Kutta"
case ("rk3.4", "RK3.4", "ssprk(4,3)", "SSPRK(4,3)")
evolve => evolve_ssprk34
order = 3
registers = 2
cfl = 2.0d+00 * cfl
name_int = "3rd order SSPRK(4,3)"
@ -256,6 +262,7 @@ module evolution
case ("rk3.5", "RK3.5", "ssprk(5,3)", "SSPRK(5,3)")
evolve => evolve_ssprk35
order = 3
registers = 2
cfl = 2.65062919143939d+00 * cfl
name_int = "3rd order SSPRK(5,3)"
@ -264,6 +271,7 @@ module evolution
error_control = .true.
evolve => evolve_ssp324
order = 3
registers = 3
stages = 4
cfl = 2.0d+00 * cfl
@ -281,6 +289,7 @@ module evolution
error_control = .true.
evolve => evolve_ssprk32m
order = 3
registers = 4
n = 2
do while(n**2 <= stages)
@ -312,6 +321,7 @@ module evolution
case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)")
evolve => evolve_ssprk4_10
order = 4
registers = 2
cfl = 6.0d+00 * cfl
name_int = "Optimal 4th order SSPRK(10,4)"
@ -591,6 +601,197 @@ module evolution
!
!===============================================================================
!
! subroutine INITIALIZE_TIME_STEP:
! -------------------------------
!
! Subroutine estimates the initial time step.
!
! References:
!
! [1] Courant, R., Friedrichs, K., Lewy, H.,
! "Über die partiellen Differenzengleichungen der mathematischen Physik",
! Mathematische Annalen, 1928, volume 100, pages 3274,
! DOI: 10.1007/BF01448839
! [2] Hairer, E., Nørsett, S. P., Wanner, G.,
! "Solving Ordinary Differential Equations I: Nonstiff Problems",
! Springer-Verlag Berlin Heidelberg, 1993,
! ISBN: 978-3-540-56670-0, DOI: 10.1007/978-3-540-78862-1
!
!===============================================================================
!
subroutine initialize_time_step()
use blocks , only : block_data, list_data, get_nleafs
use coordinates, only : nc => ncells, nb, ne, adx, ady
#if NDIMS == 3
use coordinates, only : adz
#endif /* NDIMS == 3 */
use equations , only : nf, maxspeed, cmax, cmax2
#ifdef MPI
use mpitools , only : reduce_maximum, reduce_sum
#endif /* MPI */
use sources , only : viscosity, resistivity, update_sources
implicit none
type(block_data), pointer :: pdata
integer :: mlev
real(kind=8) :: dx_min, fnorm, h0, h1
real(kind=8), dimension(3) :: d
#if NDIMS == 3
real(kind=8), dimension(nf,nc,nc,nc) :: sc, df
#else /* NDIMS == 3 */
real(kind=8), dimension(nf,nc,nc, 1) :: sc, df
#endif /* NDIMS == 3 */
real(kind=8), parameter :: eps = tiny(cmax)
!
!-------------------------------------------------------------------------------
!
#ifdef PROFILE
call start_timer(imi)
#endif /* PROFILE */
cmax = eps
mlev = 1
pdata => list_data
do while (associated(pdata))
mlev = max(mlev, pdata%meta%level)
cmax = max(cmax, maxspeed(pdata%q(:,:,:,:)))
pdata => pdata%next
end do
#ifdef MPI
call reduce_maximum(cmax)
call reduce_maximum(mlev)
#endif /* MPI */
cmax2 = cmax * cmax
#if NDIMS == 2
dx_min = min(adx(mlev), ady(mlev))
#endif /* NDIMS == 2 */
#if NDIMS == 3
dx_min = min(adx(mlev), ady(mlev), adz(mlev))
#endif /* NDIMS == 3 */
! determine the time step due to the stability
!
dth = dx_min / max(cmax &
+ 2.0d+00 * max(viscosity, resistivity) / dx_min, eps)
dt = cfl * dth
! determine the time step due to the error
!
if (error_control) then
fnorm = 1.0d+00 / (get_nleafs() * nc**NDIMS)
d(:) = 0.0d+00
pdata => list_data
do while (associated(pdata))
call update_increment(pdata)
call update_sources(pdata, 0.0d+00, 0.0d+00, pdata%du(:,:,:,:))
#if NDIMS == 3
sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))
d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne,nb:ne,1) / sc(:,:,:,:))**2)
d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne,nb:ne) / sc(:,:,:,:))**2)
#else /* NDIMS == 3 */
sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1))
d(1) = d(1) + sum((pdata%uu(:,nb:ne,nb:ne, : ,1) / sc(:,:,:,:))**2)
d(2) = d(2) + sum((pdata%du(:,nb:ne,nb:ne, : ) / sc(:,:,:,:))**2)
#endif /* NDIMS == 3 */
pdata => pdata%next
end do
#ifdef MPI
call reduce_sum(d(1:2))
#endif /* MPI */
d(1:2) = sqrt(fnorm * d(1:2))
if (minval(d(1:2)) >= 1.0d-05) then
h0 = 1.0d-02 * d(1) / d(2)
else
h0 = 1.0d-06
end if
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + h0 * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next
end do
call update_variables(time + h0, h0)
pdata => list_data
do while (associated(pdata))
#if NDIMS == 3
df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne,nb:ne)
sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))
#else /* NDIMS == 3 */
df(:,:,:,:) = pdata%du(:,nb:ne,nb:ne, : )
sc(:,:,:,:) = atol + rtol * abs(pdata%uu(:,nb:ne,nb:ne, : ,1))
#endif /* NDIMS == 3 */
call update_increment(pdata)
call update_sources(pdata, time + h0, 0.0d+00, pdata%du(:,:,:,:))
#if NDIMS == 3
df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne,nb:ne)
#else /* NDIMS == 3 */
df(:,:,:,:) = df(:,:,:,:) - pdata%du(:,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */
d(3) = d(3) + sum((df(:,:,:,:) / sc(:,:,:,:))**2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef MPI
call reduce_sum(d(3:3))
#endif /* MPI */
d(3) = sqrt(fnorm * d(3)) / h0
if (max(d(2), d(3)) > 1.0d-15) then
h1 = (1.0d-02 / max(d(2), d(3)))**(1.0d+00 / (order + 1))
else
h1 = max(1.0d-06, 1.0d-03 * h0)
end if
dte = min(1.0d+02 * h0, h1)
dt = min(dt, dte)
end if
#ifdef PROFILE
call stop_timer(imi)
#endif /* PROFILE */
!-------------------------------------------------------------------------------
!
end subroutine initialize_time_step
!
!===============================================================================
!
! subroutine NEW_TIME_STEP:
! ------------------------
!