gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2021-11-10 10:03:11 -03:00
commit cbbec5a08a
1 changed files with 581 additions and 419 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

402
sources/evolution.F90
View File

@ -230,7 +230,7 @@ module evolution
evolve => evolve_euler evolve => evolve_euler
order = 1 order = 1
registers = 1 registers = 2
name_int = "1st order Euler" name_int = "1st order Euler"
case ("rk2", "RK2") case ("rk2", "RK2")
@ -259,7 +259,7 @@ module evolution
evolve => evolve_ssprk35 evolve => evolve_ssprk35
order = 3 order = 3
registers = 2 registers = 3
cfl = 2.65062919143939d+00 * cfl cfl = 2.65062919143939d+00 * cfl
name_int = "3rd order SSPRK(5,3)" name_int = "3rd order SSPRK(5,3)"
@ -267,7 +267,7 @@ module evolution
evolve => evolve_ssprk4_10 evolve => evolve_ssprk4_10
order = 4 order = 4
registers = 2 registers = 3
cfl = 6.0d+00 * cfl cfl = 6.0d+00 * cfl
name_int = "Optimal 4th order SSPRK(10,4)" name_int = "Optimal 4th order SSPRK(10,4)"
@ -1178,6 +1178,7 @@ module evolution
use coordinates , only : adz use coordinates , only : adz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
use equations , only : nf, maxspeed, cmax, cmax2 use equations , only : nf, maxspeed, cmax, cmax2
use iso_fortran_env, only : error_unit
#ifdef MPI #ifdef MPI
use mpitools , only : reduce_maximum, reduce_sum use mpitools , only : reduce_maximum, reduce_sum
#endif /* MPI */ #endif /* MPI */
@ -1187,7 +1188,8 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: mlev, status logical :: flag
integer :: mlev, n, status
real(kind=8) :: dx_min, fnorm, h0, h1 real(kind=8) :: dx_min, fnorm, h0, h1
real(kind=8), dimension(3) :: d real(kind=8), dimension(3) :: d
@ -1198,6 +1200,8 @@ module evolution
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
real(kind=8), parameter :: eps = tiny(cmax) real(kind=8), parameter :: eps = tiny(cmax)
character(len=*), parameter :: loc = 'EVOLUTION:initialize_time_step()'
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -1278,6 +1282,11 @@ module evolution
h0 = 1.0d-06 h0 = 1.0d-06
end if end if
n = 10
flag = .true.
do while(flag)
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
@ -1290,6 +1299,26 @@ module evolution
call update_variables(time + h0, h0, status) call update_variables(time + h0, h0, status)
if (status /= 0) then
h0 = 2.5d-01 * h0
flag = n > 0
n = n - 1
else
flag = .false.
end if
end do
if (status /= 0) then
write(error_unit,"('[', a, ']: ', a)") trim(loc), &
"Could not estimate the initial time step due to " // &
"the error. Setting it to the CFL time step."
dte = dt
else
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
@ -1333,6 +1362,8 @@ module evolution
end if end if
end if
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imi) call stop_timer(imi)
#endif /* PROFILE */ #endif /* PROFILE */
@ -1463,7 +1494,7 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: status integer :: status
real(kind=8) :: tm, dtm real(kind=8) :: t
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -1471,14 +1502,27 @@ module evolution
call start_timer(imu) call start_timer(imu)
#endif /* PROFILE */ #endif /* PROFILE */
tm = time + dt status = 1
dtm = dt
do while(status /= 0)
t = time + dt
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next
end do
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, dt, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1488,9 +1532,7 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dt * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next pdata => pdata%next
end do end do
@ -1507,7 +1549,21 @@ module evolution
end do end do
end if end if
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
if (status /= 0) dt = 2.5d-01 * dt
end do
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imu) call stop_timer(imu)
@ -1545,7 +1601,7 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: status integer :: status
real(kind=8) :: tm, dtm real(kind=8) :: t, ds
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -1553,16 +1609,21 @@ module evolution
call start_timer(imu) call start_timer(imu)
#endif /* PROFILE */ #endif /* PROFILE */
status = 1
do while(status /= 0)
ds = 5.0d-01 * dt
!= 1st step: U(1) = U(n) + dt * F[U(n)] != 1st step: U(1) = U(n) + dt * F[U(n)]
! !
tm = time + dt t = time + dt
dtm = dt
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1579,18 +1640,17 @@ module evolution
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
if (status /= 0) go to 100
!= 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)] != 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)]
! !
tm = time + dt
dtm = 0.5d+00 * dt
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1600,11 +1660,9 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = 0.5d+00 * (pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = 5.0d-01 * (pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2) & + pdata%uu(:,:,:,:,2)) &
+ dt * pdata%du(:,:,:,:)) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next pdata => pdata%next
end do end do
@ -1621,7 +1679,23 @@ module evolution
end do end do
end if end if
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imu) call stop_timer(imu)
@ -1645,7 +1719,6 @@ module evolution
! "Numerical Recipes in Fortran", ! "Numerical Recipes in Fortran",
! Cambridge University Press, Cambridge, 1992 ! Cambridge University Press, Cambridge, 1992
! !
!
!=============================================================================== !===============================================================================
! !
subroutine evolve_rk3() subroutine evolve_rk3()
@ -1660,8 +1733,7 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: status integer :: status
real(kind=8) :: ds real(kind=8) :: t, ds
real(kind=8) :: tm, dtm
real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00 real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00
real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00 real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00
@ -1672,17 +1744,19 @@ module evolution
call start_timer(imu) call start_timer(imu)
#endif /* PROFILE */ #endif /* PROFILE */
status = 1
do while(status /= 0)
!= 1st substep: U(1) = U(n) + dt F[U(n)] != 1st substep: U(1) = U(n) + dt F[U(n)]
! !
ds = dt t = time + dt
tm = time + ds
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1692,26 +1766,27 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dt * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
if (status /= 0) go to 100
!= 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)] != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)]
! !
ds = f22 * dt ds = f22 * dt
tm = time + 0.5d+00 * dt t = time + 0.5d+00 * dt
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1727,19 +1802,20 @@ module evolution
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)] != 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)]
! !
ds = f32 * dt ds = f32 * dt
tm = time + dt t = time + dt
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1749,7 +1825,7 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = f31 * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,1)
@ -1769,7 +1845,23 @@ module evolution
end do end do
end if end if
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imu) call stop_timer(imu)
@ -1809,8 +1901,7 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: status integer :: status
real(kind=8) :: ds real(kind=8) :: t, ds
real(kind=8) :: tm, dtm
! !
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
@ -1818,17 +1909,20 @@ module evolution
call start_timer(imu) call start_timer(imu)
#endif /* PROFILE */ #endif /* PROFILE */
status = 1
do while(status /= 0)
!= 1st step: U(1) = U(n) + 1/2 dt F[U(n)] != 1st step: U(1) = U(n) + 1/2 dt F[U(n)]
! !
ds = dt / 2.0d+00 ds = 5.0d-01 * dt
tm = time + ds t = time + ds
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1838,24 +1932,26 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2) pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 2nd step: U(2) = U(2) + 1/2 dt F[U(1)] != 2nd step: U(2) = U(2) + 1/2 dt F[U(1)]
! !
tm = time + dt t = time + dt
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1865,22 +1961,24 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dtm * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)]) != 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)])
! !
tm = time + ds t = time + ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1892,22 +1990,24 @@ module evolution
pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2) & + pdata%uu(:,:,:,:,2) &
+ dtm * pdata%du(:,:,:,:)) / 3.0d+00 + ds * pdata%du(:,:,:,:)) / 3.0d+00
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= the final step: U(n+1) = U(3) + 1/2 dt F[U(3)] != the final step: U(n+1) = U(3) + 1/2 dt F[U(3)]
! !
tm = time + dt t = time + dt
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -1917,9 +2017,7 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next pdata => pdata%next
end do end do
@ -1936,7 +2034,23 @@ module evolution
end do end do
end if end if
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imu) call stop_timer(imu)
@ -1976,8 +2090,7 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: status integer :: status
real(kind=8) :: ds real(kind=8) :: t, ds
real(kind=8) :: tm, dtm
real(kind=8), parameter :: b1 = 3.77268915331368d-01 real(kind=8), parameter :: b1 = 3.77268915331368d-01
real(kind=8), parameter :: b3 = 2.42995220537396d-01 real(kind=8), parameter :: b3 = 2.42995220537396d-01
@ -1996,17 +2109,20 @@ module evolution
call start_timer(imu) call start_timer(imu)
#endif /* PROFILE */ #endif /* PROFILE */
status = 1
do while(status /= 0)
!= 1st step: U(1) = U(n) + b1 dt F[U(n)] != 1st step: U(1) = U(n) + b1 dt F[U(n)]
! !
ds = b1 * dt ds = b1 * dt
tm = time + ds t = time + ds
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2023,18 +2139,19 @@ module evolution
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 2nd step: U(2) = U(1) + b1 dt F[U(1)] != 2nd step: U(2) = U(1) + b1 dt F[U(1)]
! !
tm = time + 2.0d+00 * ds t = time + 2.0d+00 * ds
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2049,19 +2166,20 @@ module evolution
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] != 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)]
! !
ds = b3 * dt ds = b3 * dt
tm = time + (2.0d+00 * a33 * b1 + b3) * dt t = time + (2.0d+00 * a33 * b1 + b3) * dt
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2077,19 +2195,20 @@ module evolution
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] != 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)]
! !
ds = b4 * dt ds = b4 * dt
tm = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt t = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * dt
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2099,27 +2218,28 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = a41 * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,3) = a41 * pdata%uu(:,:,:,:,1) &
+ a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,3)
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds, status)
if (status /= 0) go to 100
!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] != the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)]
! !
ds = b5 * dt ds = b5 * dt
tm = time + dt t = time + dt
dtm = ds
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2129,8 +2249,10 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = a53 * pdata%uu(:,:,:,:,2) & pdata%uu(:,:,:,:,2) = a53 * pdata%uu(:,:,:,:,2) &
+ a55 * pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) + a55 * pdata%uu(:,:,:,:,3) + ds * pdata%du(:,:,:,:)
pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next pdata => pdata%next
end do end do
@ -2147,7 +2269,23 @@ module evolution
end do end do
end if end if
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imu) call stop_timer(imu)
@ -2185,8 +2323,8 @@ module evolution
type(block_data), pointer :: pdata type(block_data), pointer :: pdata
integer :: n, status integer :: i, status
real(kind=8) :: tm, dtm real(kind=8) :: t
real(kind=8), dimension(9) :: ds real(kind=8), dimension(9) :: ds
@ -2202,28 +2340,34 @@ module evolution
ds(:) = c(:) * dt ds(:) = c(:) * dt
status = 1
do while(status /= 0)
!= 1st step: U(2) = U(1) != 1st step: U(2) = U(1)
! !
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1)
pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1)
pdata%u => pdata%uu(:,:,:,:,2)
pdata => pdata%next pdata => pdata%next
end do end do
!= 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5 != 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5
! !
do n = 1, 5 do i = 1, 5
tm = time + ds(n) t = time + ds(i)
dtm = ds(1)
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, 0.0d+00, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2233,24 +2377,26 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:)
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds(1), status)
end do ! n = 1, 5 if (status /= 0) go to 100
end do ! i = 1, 5
!= 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1) != 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1)
! !
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,2) = (pdata%uu(:,:,:,:,2) & pdata%uu(:,:,:,:,3) = (pdata%uu(:,:,:,:,3) &
+ 9.0d+00 * pdata%uu(:,:,:,:,1)) / 2.5d+01 + 9.0d+00 * pdata%uu(:,:,:,:,2)) / 2.5d+01
pdata%uu(:,:,:,:,1) = 1.5d+01 * pdata%uu(:,:,:,:,2) & pdata%uu(:,:,:,:,2) = 1.5d+01 * pdata%uu(:,:,:,:,3) &
- 5.0d+00 * pdata%uu(:,:,:,:,1) - 5.0d+00 * pdata%uu(:,:,:,:,2)
pdata => pdata%next pdata => pdata%next
end do end do
@ -2259,16 +2405,15 @@ module evolution
! !
! integrate the intermediate steps ! integrate the intermediate steps
! !
do n = 6, 9 do i = 6, 9
tm = time + ds(n) t = time + ds(i)
dtm = ds(1)
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, ds(1), pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2278,25 +2423,26 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds(1) * pdata%du(:,:,:,:)
pdata => pdata%next pdata => pdata%next
end do end do
call update_variables(tm, dtm, status) call update_variables(t, ds(1), status)
end do ! n = 6, 9 if (status /= 0) go to 100
end do ! i = 6, 9
!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] != the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)]
! !
tm = time + dt t = time + dt
dtm = dt / 1.0d+01
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
call update_increment(pdata) call update_increment(pdata)
call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) call update_sources(pdata, t, 1.0d-01 * dt, pdata%du(:,:,:,:))
pdata => pdata%next pdata => pdata%next
end do end do
@ -2306,9 +2452,9 @@ module evolution
pdata => list_data pdata => list_data
do while (associated(pdata)) do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) & pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,3) &
+ 6.0d-01 * pdata%uu(:,:,:,:,1) & + 6.0d-01 * pdata%uu(:,:,:,:,2) &
+ dtm * pdata%du(:,:,:,:) + 1.0d-01 * dt * pdata%du(:,:,:,:)
pdata => pdata%next pdata => pdata%next
end do end do
@ -2325,7 +2471,23 @@ module evolution
end do end do
end if end if
call update_variables(tm, dtm, status) call update_variables(t, dt, status)
100 continue
if (status /= 0) dt = 2.5d-01 * dt
end do
pdata => list_data
do while (associated(pdata))
pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2)
pdata%u => pdata%uu(:,:,:,:,1)
pdata => pdata%next
end do
#ifdef PROFILE #ifdef PROFILE
call stop_timer(imu) call stop_timer(imu)