EVOLUTION: Rewrite SSPRK4(2)M integration method.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
e7612f5cfd
commit
0c8d8fac9c
@ -276,11 +276,11 @@ module evolution
|
||||
|
||||
betas(:) = - betas(:) / 3.0d+00
|
||||
|
||||
case ("rk3.m", "ssprk(m,3)", "SSPRK(m,3)")
|
||||
case ("SSPRK3(2)m", "SSP3(2)m", "ssprk3(2)m", "ssp3(2)m", "SSPRK(m,3)", "ssprk(m,3)")
|
||||
|
||||
error_control = .true.
|
||||
evolve => evolve_ssprk3_m
|
||||
registers = 3
|
||||
evolve => evolve_ssprk32m
|
||||
registers = 4
|
||||
n = 2
|
||||
do while(n**2 <= stages)
|
||||
n = n + 1
|
||||
@ -288,7 +288,25 @@ module evolution
|
||||
n = n - 1
|
||||
stages = max(4, n**2)
|
||||
cfl = (n - 1) * n * cfl
|
||||
write(name_int, "('Optimal 3rd order SSPRK(',i0,',3) with error control')") stages
|
||||
write(name_int, "('3ʳᵈ-order order ',i0,'-step embedded SSPRK3(2)n² method')") stages
|
||||
|
||||
betas(:) = [ 5.5d-01, -2.7d-01, 5.0d-02 ]
|
||||
|
||||
call get_parameter("beta(1)", betas(1))
|
||||
call get_parameter("beta(2)", betas(2))
|
||||
call get_parameter("beta(3)", betas(3))
|
||||
|
||||
betas(:) = - betas(:) / 3.0d+00
|
||||
|
||||
if (stages == 4 .and. verbose) then
|
||||
write(*,*)
|
||||
write(*,*) "ADVICE:"
|
||||
write(*,*) "The method 'SSPRK3(2)m' is primarily designed for n²," // &
|
||||
" where n > 2. Since you intend to use 4 stages with" // &
|
||||
" this method, it is highly advised to use method" // &
|
||||
" 'SSPRK3(2)4', since it slightly better optimized and" // &
|
||||
" significantly used less memory."
|
||||
end if
|
||||
|
||||
case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)")
|
||||
|
||||
@ -307,6 +325,8 @@ module evolution
|
||||
write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," // &
|
||||
" 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," // &
|
||||
" 'ssprk(m,3)', 'ssprk(10,4)'."
|
||||
write(*,"(1x,a)") "Available embedded methods: 'ssprk3(2)4'," // &
|
||||
" 'ssprk3(2)m'."
|
||||
end if
|
||||
status = 1
|
||||
|
||||
@ -1943,30 +1963,31 @@ module evolution
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine EVOLVE_SSPRK3_M:
|
||||
! subroutine EVOLVE_SSPRK32M:
|
||||
! --------------------------
|
||||
!
|
||||
! Subroutine advances the solution by one time step using the 3rd order
|
||||
! m-stage Strong Stability Preserving Runge-Kutta time integration method.
|
||||
! Subroutine advances the solution by one time step using the 3ʳᵈ-order
|
||||
! m-stage (m=n²) embedded Strong Stability Preserving Runge-Kutta time
|
||||
! integration method SSP3(2)m with the error control.
|
||||
!
|
||||
! References:
|
||||
!
|
||||
! [1] Gottlieb, S., Ketcheson, D., Shu C. - W.,
|
||||
! "Strong stability preserving Runge-Kutta and multistep
|
||||
! time discretizations",
|
||||
! World Scientific Publishing, 2011
|
||||
! [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., Ketcheson, D., Shu, C.-W.,
|
||||
! "Strong stability preserving Runge-Kutta and multistep time
|
||||
! discretizations",
|
||||
! 2011, World Scientific Publishing
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine evolve_ssprk3_m()
|
||||
subroutine evolve_ssprk32m()
|
||||
|
||||
use blocks , only : block_data, list_data, get_dblocks
|
||||
use blocks , only : block_data, list_data
|
||||
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 equations , only : errors, ibp, cmax
|
||||
use sources , only : update_sources
|
||||
|
||||
implicit none
|
||||
@ -1974,18 +1995,12 @@ module evolution
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
logical , save :: first = .true.
|
||||
integer , save :: n, n1, n2, n3, n4
|
||||
real(kind=8), save :: r, f1, f2, g1, g2
|
||||
integer , save :: i1, i2, i3, i4, n, m
|
||||
real(kind=8), save :: f1, f2, f3, f4
|
||||
|
||||
logical :: test
|
||||
integer :: i, l, nrej
|
||||
real(kind=8) :: tm, dtm, ds, umax, emax
|
||||
real(kind=8) :: fc, fcmn, fcmx
|
||||
|
||||
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
|
||||
|
||||
real(kind=8), parameter :: k1 = -5.8d-01 / 3.0d+00, k2 = 7.0d-02, &
|
||||
k3 = -1.0d-01 / 3.0d+00
|
||||
integer :: nrej, i
|
||||
real(kind=8) :: tm, dtm, dh, fc
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -1995,45 +2010,45 @@ module evolution
|
||||
|
||||
if (first) then
|
||||
|
||||
n = int(sqrt(1.0d+00 * stages))
|
||||
n1 = (n - 1) * (n - 2) / 2
|
||||
n2 = n1 + 1
|
||||
n3 = n * (n + 1) / 2 - 1
|
||||
n4 = n * (n + 1) / 2 + 1
|
||||
n = 2
|
||||
do while(n**2 < stages)
|
||||
n = n + 1
|
||||
end do
|
||||
m = stages - n
|
||||
i1 = (n - 1) * (n - 2) / 2
|
||||
i2 = i1 + 1
|
||||
i3 = n * (n + 1) / 2
|
||||
i4 = i3 + 1
|
||||
|
||||
r = (1.0d+00 * (2 * n - 1))
|
||||
f1 = (1.0d+00 * n ) / r
|
||||
f2 = (1.0d+00 * (n - 1)) / r
|
||||
r = 1.0d+00 * (stages - n)
|
||||
|
||||
g1 = 1.0d+00 / ( stages - n) - 1.0d+00 / stages
|
||||
g2 = 1.0d+00 / (2 * stages - n) - 1.0d+00 / stages
|
||||
fc = real(2 * n - 1, kind=8)
|
||||
f1 = real(n - 1, kind=8) / fc
|
||||
f2 = real(n , kind=8) / fc
|
||||
f4 = real(m , kind=8) / stages
|
||||
f3 = 1.0d+00 - f4
|
||||
|
||||
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 3rd order solution, U(1)
|
||||
! pdata%uu(:,:,:,:,3) - the intermediate solution, U(2)
|
||||
! pdata%uu(:,:,:,:,4) - the new 2nd order solution, U(3)
|
||||
!
|
||||
do while(test)
|
||||
|
||||
lerr(:,:,:,:,:) = 0.0d+00
|
||||
! initiate the fractional time step
|
||||
!
|
||||
dh = dt / m
|
||||
|
||||
ds = dt / r
|
||||
|
||||
!= 1st step: U(1) = U(n)
|
||||
!= preparation step: U(1) = U(n)
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2045,12 +2060,12 @@ module evolution
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2
|
||||
!= 1st step: U(1) = U(1) + dt/r F[U(1)], for i = 1...(n-1)*(n-2)/2
|
||||
!
|
||||
do i = 1, n1
|
||||
do i = 1, i1
|
||||
|
||||
tm = time + i * ds
|
||||
dtm = ds
|
||||
tm = time + (i - 1) * dh
|
||||
dtm = dh
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2063,29 +2078,19 @@ 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,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! n = 1, n1
|
||||
end do ! i = 1, i1
|
||||
|
||||
!= 3rd step: U(2) = U(1)
|
||||
!
|
||||
! iterate over all data blocks
|
||||
!= 2nd step: U(2) = U(1)
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2093,15 +2098,14 @@ module evolution
|
||||
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,2)
|
||||
|
||||
pdata => pdata%next
|
||||
|
||||
end do
|
||||
|
||||
!= 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1
|
||||
!= 3rd step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1...n*(n+1)/2
|
||||
!
|
||||
do i = n2, n3
|
||||
do i = i2, i3
|
||||
|
||||
tm = time + i * ds
|
||||
dtm = ds
|
||||
tm = time + (i - 1) * dh
|
||||
dtm = dh
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2114,30 +2118,43 @@ 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,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! i = n2, n3
|
||||
end do ! i = i2, i3
|
||||
|
||||
!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1)
|
||||
!= 4th step: U(3) = (1 - r/s) * U(n) + r/s * U(1),
|
||||
! U(1) = [(n - 1) * U(1) + n * U(2)] / (2*n - 1),
|
||||
! U(3) = U(3) - r/s * U(1)
|
||||
!
|
||||
tm = time + dt
|
||||
dtm = ds
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%uu(:,:,:,:,4) = f3 * pdata%uu(:,:,:,:,1) &
|
||||
+ f4 * pdata%uu(:,:,:,:,2)
|
||||
pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,2) &
|
||||
+ f2 * pdata%uu(:,:,:,:,3)
|
||||
pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) &
|
||||
- f4 * pdata%uu(:,:,:,:,2)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2+1...stages
|
||||
!
|
||||
do i = i4, stages
|
||||
|
||||
tm = time + (i - 1 - n) * dh
|
||||
dtm = dh
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
@ -2150,131 +2167,90 @@ module evolution
|
||||
|
||||
call boundary_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,3) &
|
||||
+ f2 * (pdata%uu(:,:,:,:,2) &
|
||||
+ ds * pdata%du(:,:,:,:))
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m
|
||||
end do ! i = i4, stages
|
||||
|
||||
!= 7th step: U(3) = U(3) + r/s * U(1)
|
||||
!
|
||||
do i = n4, stages
|
||||
|
||||
tm = time + (i - n) * ds
|
||||
dtm = ds
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
call update_increment(pdata)
|
||||
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:))
|
||||
pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) + f4 * pdata%uu(:,:,:,:,2)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call boundary_fluxes()
|
||||
!= 8th step: decay the magnetic divergence potential of both solutions
|
||||
!
|
||||
if (ibp > 0) then
|
||||
adecay(:) = exp(aglm(:) * cmax * dt)
|
||||
|
||||
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,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
|
||||
* pdata%uu(ibp,:,:,:,2)
|
||||
pdata%uu(ibp,:,:,:,4) = adecay(pdata%meta%level) &
|
||||
* pdata%uu(ibp,:,:,:,4)
|
||||
|
||||
pdata => pdata%next
|
||||
end do
|
||||
end if
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! i = n4, stages
|
||||
|
||||
! find umax
|
||||
! update the vector of errors
|
||||
!
|
||||
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
|
||||
call update_errors(2, 4)
|
||||
|
||||
! get the maximum error of each variable
|
||||
! calculate the tolerance and estimate the next time step due to the error
|
||||
!
|
||||
do l = 1, nf
|
||||
errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
|
||||
end do
|
||||
errtol = maxval(errors)
|
||||
errs(1) = errtol
|
||||
fc = product(errs(:)**betas(:))
|
||||
fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi)
|
||||
dte = dt * fc
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(umax)
|
||||
call reduce_maximum(errors)
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate tolerance and time step
|
||||
!
|
||||
emax = maxval(errors) / (atol + rtol * umax)
|
||||
|
||||
if (emax <= 1.0d+00 .or. nrej >= mrej) then
|
||||
if (fc > fac .or. nrej >= mrej) 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
|
||||
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)
|
||||
|
||||
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 = dt / r
|
||||
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
@ -2285,27 +2261,13 @@ module evolution
|
||||
pdata => pdata%next
|
||||
end do
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
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
|
||||
|
||||
#ifdef PROFILE
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk3_m
|
||||
end subroutine evolve_ssprk32m
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user