EVOLUTION: Rewrite SSPRK2(1)m method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
a87bcce907
commit
a9cf789833
@ -279,7 +279,30 @@ module evolution
|
||||
registers = 3
|
||||
stages = max(2, min(9, stages))
|
||||
cfl = (stages - 1) * cfl
|
||||
write(name_int, "('Optimal 2nd order SSPRK(',i0,',2) with error control')") stages
|
||||
write(name_int, "('2ⁿᵈ-order ',i0,'-step embedded " // &
|
||||
"SSPRK(',i0,',2) method')") stages, stages
|
||||
|
||||
allocate(c(stages), dl(stages), bt(stages), bh(stages))
|
||||
|
||||
do n = 1, stages
|
||||
c(n) = n - 1.0d+00
|
||||
end do
|
||||
c (:) = c(:) / (stages - 1)
|
||||
bt(:) = 1.0d+00 / (stages - 1)
|
||||
bt(stages) = 1.0d+00 / stages
|
||||
bh(:) = 1.0d+00 / stages
|
||||
bh(1) = (stages + 1.0d+00) / stages**2
|
||||
bh(stages) = (stages - 1.0d+00) / stages**2
|
||||
dl(1) = (stages - 1.0d+00) / stages
|
||||
dl(2) = 1.0d+00 - dl(1)
|
||||
|
||||
betas(:) = [ 5.8d-01, -2.1d-01, 1.0d-01 ]
|
||||
|
||||
call get_parameter("beta(1)", betas(1))
|
||||
call get_parameter("beta(2)", betas(2))
|
||||
call get_parameter("beta(3)", betas(3))
|
||||
|
||||
betas(:) = - betas(:) / 2.0d+00
|
||||
|
||||
case ("SSPRK3(2)4", "SSP3(2)4", "ssprk3(2)4", "ssp3(2)4")
|
||||
|
||||
@ -2327,7 +2350,11 @@ module evolution
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Gottlieb, S. and Gottlieb, L.-A., J.
|
||||
! [1] Conde, S., Fekete, I., Shadid, J. N.,
|
||||
! "Embedded error estimation and adaptive step-size control for
|
||||
! optimal explicit strong stability preserving Runge–Kutta methods"
|
||||
! 2018, arXiv:1806.08693v1
|
||||
! [2] Gottlieb, S. and Gottlieb, L.-A., J.
|
||||
! "Strong Stability Preserving Properties of Runge-Kutta Time
|
||||
! Discretization Methods for Linear Constant Coefficient Operators",
|
||||
! Journal of Scientific Computing,
|
||||
@ -2337,71 +2364,43 @@ module evolution
|
||||
!
|
||||
subroutine evolve_ssprk2_m()
|
||||
|
||||
use blocks , only : block_data, list_data, get_dblocks
|
||||
use boundaries , only : boundary_fluxes
|
||||
use coordinates, only : nc => ncells, nb, ne
|
||||
use equations , only : errors, nf, ibp, cmax
|
||||
#ifdef MPI
|
||||
use mpitools , only : reduce_maximum
|
||||
#endif /* MPI */
|
||||
use sources , only : update_sources
|
||||
use blocks , only : block_data, list_data
|
||||
use boundaries , only : boundary_fluxes
|
||||
use coordinates , only : nb, ne
|
||||
use equations , only : errors, ibp, cmax
|
||||
use iso_fortran_env, only : error_unit
|
||||
use sources , only : update_sources
|
||||
|
||||
implicit none
|
||||
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
logical :: test
|
||||
integer :: n, l, nrej, status
|
||||
real(kind=8) :: tm, dtm, ds, umax, emax
|
||||
real(kind=8) :: fc, fcmn, fcmx
|
||||
logical :: test
|
||||
integer :: nrej, i, status
|
||||
real(kind=8) :: tm, dtm, dh, fc
|
||||
|
||||
logical , save :: first = .true.
|
||||
real(kind=8), save :: ft, fl, fr, gt, gl
|
||||
character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk2_m()'
|
||||
|
||||
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
|
||||
|
||||
real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
if (first) then
|
||||
|
||||
ft = 1.0d+00 / (stages - 1)
|
||||
fl = 1.0d+00 / stages
|
||||
fr = 1.0d+00 - fl
|
||||
|
||||
gt = fl - ft
|
||||
gl = fl
|
||||
|
||||
first = .false.
|
||||
|
||||
end if
|
||||
|
||||
fc = fac
|
||||
fcmn = facmin
|
||||
fcmx = facmax
|
||||
|
||||
l = get_dblocks()
|
||||
#if NDIMS == 3
|
||||
allocate(lerr(l,nf,nc,nc,nc))
|
||||
#else /* NDIMS == 3 */
|
||||
allocate(lerr(l,nf,nc,nc, 1))
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
test = .true.
|
||||
nrej = 0
|
||||
|
||||
! at the entry point we assume the previous solution of conserved variables U(n)
|
||||
! is stored in pdata%u(:,:,:,:,1) and the primitive variables are already
|
||||
! updated from this array;
|
||||
!
|
||||
! pdata%uu(:,:,:,:,1) - the previous solution, U(n)
|
||||
! pdata%uu(:,:,:,:,2) - the new 2nd order solution, U(1)
|
||||
! pdata%uu(:,:,:,:,3) - the new 1st order solution, U(2)
|
||||
!
|
||||
do while(test)
|
||||
|
||||
lerr(:,:,:,:,:) = 0.0d+00
|
||||
|
||||
ds = ft * dt
|
||||
|
||||
!= 1st step: U(1) = U(n), U(2) = U(n)
|
||||
!= preparatory step: U(1) = U(n), U(2) = 0, U(3) = U(n)
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2414,12 +2413,12 @@ module evolution
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1
|
||||
!= 1st step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1
|
||||
!
|
||||
do n = 1, stages - 1
|
||||
do i = 1, stages - 1
|
||||
|
||||
tm = time + n * ds
|
||||
dtm = ds
|
||||
dtm = c(i) * dt
|
||||
tm = time + dtm
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2432,30 +2431,27 @@ module evolution
|
||||
|
||||
call boundary_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) &
|
||||
+ bt(i) * dt * pdata%du(:,:,:,:)
|
||||
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) &
|
||||
+ bh(i) * dt * pdata%du(:,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm, status)
|
||||
if (status /= 0) go to 100
|
||||
|
||||
end do ! n = 1, stages - 1
|
||||
end do ! i = 1, stages - 1
|
||||
|
||||
!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)]
|
||||
!= 2nd step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)]
|
||||
!
|
||||
tm = time + dt
|
||||
dtm = fl * dt
|
||||
i = stages
|
||||
dtm = c(i) * dt
|
||||
tm = time + dtm
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2468,93 +2464,93 @@ module evolution
|
||||
|
||||
call boundary_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%uu(:,:,:,:,2) = fr * pdata%uu(:,:,:,:,2) &
|
||||
+ fl * (pdata%uu(:,:,:,:,3) &
|
||||
+ dt * pdata%du(:,:,:,:))
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
pdata%uu(:,:,:,:,2) = dl(1) * pdata%uu(:,:,:,:,2) &
|
||||
+ dl(2) * pdata%uu(:,:,:,:,1) &
|
||||
+ bt(i) * dt * pdata%du(:,:,:,:)
|
||||
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) &
|
||||
+ bh(i) * dt * pdata%du(:,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm, status)
|
||||
if (status /= 0) go to 100
|
||||
|
||||
! find umax
|
||||
!= 3rd step: decay the magnetic divergence potential of both solutions
|
||||
!
|
||||
umax = 0.0d+00
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
#if NDIMS == 3
|
||||
umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))), &
|
||||
maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,2))))
|
||||
#else /* NDIMS == 3 */
|
||||
umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,1))), &
|
||||
maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,2))))
|
||||
#endif /* NDIMS == 3 */
|
||||
pdata => pdata%next
|
||||
end do
|
||||
if (ibp > 0) then
|
||||
adecay(:) = exp(aglm(:) * cmax * dt)
|
||||
|
||||
! get the maximum error of each variable
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
|
||||
* pdata%uu(ibp,:,:,:,2)
|
||||
pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level) &
|
||||
* pdata%uu(ibp,:,:,:,3)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
end if
|
||||
|
||||
! update the vector of errors
|
||||
!
|
||||
do l = 1, nf
|
||||
errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
|
||||
end do
|
||||
call update_errors(2, 3)
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(umax)
|
||||
call reduce_maximum(errors)
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate tolerance and time step
|
||||
! calculate the tolerance and estimate the next time step due to the error
|
||||
!
|
||||
emax = maxval(errors) / (atol + rtol * umax)
|
||||
errtol = maxval(errors)
|
||||
errs(1) = errtol
|
||||
fc = product(errs(:)**betas(:))
|
||||
fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi)
|
||||
dte = dt * fc
|
||||
|
||||
if (emax <= 1.0d+00 .or. nrej >= mrej) then
|
||||
100 continue
|
||||
|
||||
if ((fc > fac .or. nrej >= mrej) .and. status == 0) then
|
||||
test = .false.
|
||||
|
||||
errs(3) = errs(2)
|
||||
errs(2) = errs(1)
|
||||
errs(1) = emax
|
||||
|
||||
errtol = emax
|
||||
|
||||
dte = dt * min(fcmx, max(fcmn, &
|
||||
fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
|
||||
|
||||
fcmx = facmax
|
||||
else
|
||||
errs(1) = emax
|
||||
|
||||
dte = dt * min(fcmx, max(fcmn, &
|
||||
fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
|
||||
dt = dte
|
||||
|
||||
fcmx = fac
|
||||
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
if (status == 0) then
|
||||
dt = dte
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
else
|
||||
dt = 2.5d-01 * dt
|
||||
end if
|
||||
nrejections = nrejections + 1
|
||||
|
||||
! since the solution was rejected, we have to revert the primitive variables
|
||||
! to the previous state
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%u => pdata%uu(:,:,:,:,1)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm, status)
|
||||
if (status /= 0) then
|
||||
write(error_unit,"('[', a, ']: ', a)") trim(loc), &
|
||||
"Possible bug: the revert to the previous state " //&
|
||||
"found it unphysical."
|
||||
end if
|
||||
|
||||
end if
|
||||
|
||||
niterations = niterations + 1
|
||||
|
||||
end do
|
||||
|
||||
if (allocated(lerr)) deallocate(lerr)
|
||||
|
||||
!= final step: U(n+1) = U(1)
|
||||
!= final step: U(n+1) = U(1) - update the accepted solution
|
||||
!
|
||||
tm = time + dt
|
||||
dtm = ft * dt
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
@ -2565,20 +2561,6 @@ module evolution
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
if (ibp > 0) then
|
||||
adecay(:) = exp(aglm(:) * cmax * dt)
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
end if
|
||||
|
||||
call update_variables(tm, dtm, status)
|
||||
|
||||
#ifdef PROFILE
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user