Merge branch 'master' into reconnection
This commit is contained in:
commit
090d6a6896
@ -279,7 +279,30 @@ module evolution
|
|||||||
registers = 3
|
registers = 3
|
||||||
stages = max(2, min(9, stages))
|
stages = max(2, min(9, stages))
|
||||||
cfl = (stages - 1) * cfl
|
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")
|
case ("SSPRK3(2)4", "SSP3(2)4", "ssprk3(2)4", "ssp3(2)4")
|
||||||
|
|
||||||
@ -2327,7 +2350,11 @@ module evolution
|
|||||||
!
|
!
|
||||||
! References:
|
! 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
|
! "Strong Stability Preserving Properties of Runge-Kutta Time
|
||||||
! Discretization Methods for Linear Constant Coefficient Operators",
|
! Discretization Methods for Linear Constant Coefficient Operators",
|
||||||
! Journal of Scientific Computing,
|
! Journal of Scientific Computing,
|
||||||
@ -2337,71 +2364,43 @@ module evolution
|
|||||||
!
|
!
|
||||||
subroutine evolve_ssprk2_m()
|
subroutine evolve_ssprk2_m()
|
||||||
|
|
||||||
use blocks , only : block_data, list_data, get_dblocks
|
use blocks , only : block_data, list_data
|
||||||
use boundaries , only : boundary_fluxes
|
use boundaries , only : boundary_fluxes
|
||||||
use coordinates, only : nc => ncells, nb, ne
|
use coordinates , only : nb, ne
|
||||||
use equations , only : errors, nf, ibp, cmax
|
use equations , only : errors, ibp, cmax
|
||||||
#ifdef MPI
|
use iso_fortran_env, only : error_unit
|
||||||
use mpitools , only : reduce_maximum
|
use sources , only : update_sources
|
||||||
#endif /* MPI */
|
|
||||||
use sources , only : update_sources
|
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
type(block_data), pointer :: pdata
|
type(block_data), pointer :: pdata
|
||||||
|
|
||||||
logical :: test
|
logical :: test
|
||||||
integer :: n, l, nrej, status
|
integer :: nrej, i, status
|
||||||
real(kind=8) :: tm, dtm, ds, umax, emax
|
real(kind=8) :: tm, dtm, dh, fc
|
||||||
real(kind=8) :: fc, fcmn, fcmx
|
|
||||||
|
|
||||||
logical , save :: first = .true.
|
character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk2_m()'
|
||||||
real(kind=8), save :: ft, fl, fr, gt, gl
|
|
||||||
|
|
||||||
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
|
|
||||||
|
|
||||||
real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02
|
|
||||||
!
|
|
||||||
!-------------------------------------------------------------------------------
|
!-------------------------------------------------------------------------------
|
||||||
!
|
!
|
||||||
#ifdef PROFILE
|
#ifdef PROFILE
|
||||||
call start_timer(imu)
|
call start_timer(imu)
|
||||||
#endif /* PROFILE */
|
#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.
|
test = .true.
|
||||||
nrej = 0
|
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)
|
do while(test)
|
||||||
|
|
||||||
lerr(:,:,:,:,:) = 0.0d+00
|
!= preparatory step: U(1) = U(n), U(2) = 0, U(3) = U(n)
|
||||||
|
|
||||||
ds = ft * dt
|
|
||||||
|
|
||||||
!= 1st step: U(1) = U(n), U(2) = U(n)
|
|
||||||
!
|
!
|
||||||
pdata => list_data
|
pdata => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pdata))
|
||||||
@ -2414,12 +2413,12 @@ module evolution
|
|||||||
pdata => pdata%next
|
pdata => pdata%next
|
||||||
end do
|
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 = c(i) * dt
|
||||||
dtm = ds
|
tm = time + dtm
|
||||||
|
|
||||||
pdata => list_data
|
pdata => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pdata))
|
||||||
@ -2432,30 +2431,27 @@ module evolution
|
|||||||
|
|
||||||
call boundary_fluxes()
|
call boundary_fluxes()
|
||||||
|
|
||||||
l = 1
|
|
||||||
pdata => list_data
|
pdata => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pdata))
|
||||||
|
|
||||||
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
|
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) &
|
||||||
|
+ bt(i) * dt * pdata%du(:,:,:,:)
|
||||||
#if NDIMS == 3
|
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) &
|
||||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne,nb:ne)
|
+ bh(i) * dt * pdata%du(:,:,:,:)
|
||||||
#else /* NDIMS == 3 */
|
|
||||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne, : )
|
|
||||||
#endif /* NDIMS == 3 */
|
|
||||||
l = l + 1
|
|
||||||
|
|
||||||
pdata => pdata%next
|
pdata => pdata%next
|
||||||
end do
|
end do
|
||||||
|
|
||||||
call update_variables(tm, dtm, status)
|
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
|
i = stages
|
||||||
dtm = fl * dt
|
dtm = c(i) * dt
|
||||||
|
tm = time + dtm
|
||||||
|
|
||||||
pdata => list_data
|
pdata => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pdata))
|
||||||
@ -2468,93 +2464,93 @@ module evolution
|
|||||||
|
|
||||||
call boundary_fluxes()
|
call boundary_fluxes()
|
||||||
|
|
||||||
l = 1
|
|
||||||
pdata => list_data
|
pdata => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pdata))
|
||||||
|
|
||||||
pdata%uu(:,:,:,:,2) = fr * pdata%uu(:,:,:,:,2) &
|
pdata%uu(:,:,:,:,2) = dl(1) * pdata%uu(:,:,:,:,2) &
|
||||||
+ fl * (pdata%uu(:,:,:,:,3) &
|
+ dl(2) * pdata%uu(:,:,:,:,1) &
|
||||||
+ dt * pdata%du(:,:,:,:))
|
+ bt(i) * dt * pdata%du(:,:,:,:)
|
||||||
|
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) &
|
||||||
#if NDIMS == 3
|
+ bh(i) * dt * pdata%du(:,:,:,:)
|
||||||
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 => pdata%next
|
pdata => pdata%next
|
||||||
end do
|
end do
|
||||||
|
|
||||||
call update_variables(tm, dtm, status)
|
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
|
if (ibp > 0) then
|
||||||
pdata => list_data
|
adecay(:) = exp(aglm(:) * cmax * dt)
|
||||||
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
|
|
||||||
|
|
||||||
! 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
|
call update_errors(2, 3)
|
||||||
errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
|
|
||||||
end do
|
|
||||||
|
|
||||||
#ifdef MPI
|
! calculate the tolerance and estimate the next time step due to the error
|
||||||
call reduce_maximum(umax)
|
|
||||||
call reduce_maximum(errors)
|
|
||||||
#endif /* MPI */
|
|
||||||
|
|
||||||
! calculate tolerance and time step
|
|
||||||
!
|
!
|
||||||
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.
|
test = .false.
|
||||||
|
|
||||||
errs(3) = errs(2)
|
errs(3) = errs(2)
|
||||||
errs(2) = errs(1)
|
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
|
else
|
||||||
errs(1) = emax
|
if (status == 0) then
|
||||||
|
dt = dte
|
||||||
dte = dt * min(fcmx, max(fcmn, &
|
nrej = nrej + 1 ! rejection count in the current step
|
||||||
fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3))
|
else
|
||||||
dt = dte
|
dt = 2.5d-01 * dt
|
||||||
|
end if
|
||||||
fcmx = fac
|
|
||||||
|
|
||||||
nrej = nrej + 1 ! rejection count in the current step
|
|
||||||
nrejections = nrejections + 1
|
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
|
end if
|
||||||
|
|
||||||
niterations = niterations + 1
|
niterations = niterations + 1
|
||||||
|
|
||||||
end do
|
end do
|
||||||
|
|
||||||
if (allocated(lerr)) deallocate(lerr)
|
!= final step: U(n+1) = U(1) - update the accepted solution
|
||||||
|
|
||||||
!= final step: U(n+1) = U(1)
|
|
||||||
!
|
!
|
||||||
tm = time + dt
|
|
||||||
dtm = ft * dt
|
|
||||||
|
|
||||||
pdata => list_data
|
pdata => list_data
|
||||||
do while (associated(pdata))
|
do while (associated(pdata))
|
||||||
|
|
||||||
@ -2565,20 +2561,6 @@ module evolution
|
|||||||
pdata => pdata%next
|
pdata => pdata%next
|
||||||
end do
|
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
|
#ifdef PROFILE
|
||||||
call stop_timer(imu)
|
call stop_timer(imu)
|
||||||
#endif /* PROFILE */
|
#endif /* PROFILE */
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user