Merge branch 'master' into reconnection
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
commit
7b968dc728
@ -260,10 +260,11 @@ module blocks
|
||||
real(kind=8), dimension(:,:,:,:) , pointer :: u
|
||||
|
||||
! an allocatable arrays to store all conserved
|
||||
! variables (required two for Runge-Kutta
|
||||
! temporal integration methods)
|
||||
! variables (requires two additional arrays
|
||||
! for embedded low-storage explicit optimal
|
||||
! high-order SSPRK integration methods)
|
||||
!
|
||||
real(kind=8), dimension(:,:,:,:) , allocatable :: u0, u1
|
||||
real(kind=8), dimension(:,:,:,:) , allocatable :: u0, u1, u2
|
||||
|
||||
! an allocatable array to store all primitive
|
||||
! variables
|
||||
@ -1178,8 +1179,8 @@ module blocks
|
||||
! allocate space for conserved and primitive variables, and fluxes
|
||||
!
|
||||
allocate(pdata%u0(nvars,nx,ny,nz), pdata%u1(nvars,nx,ny,nz), &
|
||||
pdata%q (nvars,nx,ny,nz), pdata%f(ndims,nflux,nx,ny,nz), &
|
||||
stat = status)
|
||||
pdata%u2(nvars,nx,ny,nz), pdata%q (nvars,nx,ny,nz), &
|
||||
pdata%f(ndims,nflux,nx,ny,nz), stat = status)
|
||||
|
||||
! check if allocation succeeded
|
||||
!
|
||||
@ -1189,6 +1190,7 @@ module blocks
|
||||
!
|
||||
pdata%u0 = 0.0d+00
|
||||
pdata%u1 = 0.0d+00
|
||||
pdata%u2 = 0.0d+00
|
||||
pdata%q = 0.0d+00
|
||||
pdata%f = 0.0d+00
|
||||
|
||||
@ -1282,6 +1284,7 @@ module blocks
|
||||
!
|
||||
if (allocated(pdata%u0)) deallocate(pdata%u0, stat = status)
|
||||
if (allocated(pdata%u1)) deallocate(pdata%u1, stat = status)
|
||||
if (allocated(pdata%u2)) deallocate(pdata%u2, stat = status)
|
||||
|
||||
! deallocate primitive variables
|
||||
!
|
||||
|
||||
@ -5805,42 +5805,47 @@ module boundaries
|
||||
!
|
||||
do p = 1, nv
|
||||
|
||||
dq(1) = 1.25d-01 * (qn(p,i+1,j,k) - qn(p,i-1,j,k))
|
||||
dq(2) = 1.25d-01 * (qn(p,i,j+1,k) - qn(p,i,j-1,k))
|
||||
dq(3) = 1.25d-01 * (qn(p,i,j,k+1) - qn(p,i,j,k-1))
|
||||
|
||||
if (positive(p) .and. qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) then
|
||||
|
||||
if (qn(p,i,j,k) <= 0.0d+00) then
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
|
||||
dq(:) = 0.0d+00
|
||||
else
|
||||
|
||||
! calculate limited derivatives in all directions
|
||||
!
|
||||
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
|
||||
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
|
||||
dq(1) = limiter_prol(0.5d+00, dql, dqr)
|
||||
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
|
||||
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
|
||||
dq(1) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
|
||||
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
|
||||
dq(2) = limiter_prol(0.5d+00, dql, dqr)
|
||||
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
|
||||
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
|
||||
dq(2) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
|
||||
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
|
||||
dq(3) = limiter_prol(0.5d+00, dql, dqr)
|
||||
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
|
||||
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
|
||||
dq(3) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
if (positive(p)) then
|
||||
if (qn(p,i,j,k) > 0.0d+00) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
|
||||
! calculate the derivative terms
|
||||
! calculate derivative terms
|
||||
!
|
||||
dq1 = dq(1) + dq(2) + dq(3)
|
||||
dq2 = dq(1) - dq(2) - dq(3)
|
||||
dq3 = dq(1) - dq(2) + dq(3)
|
||||
dq4 = dq(1) + dq(2) - dq(3)
|
||||
|
||||
! prolong the face region to the output array
|
||||
! prolong the face region
|
||||
!
|
||||
qb(p,is,js,ks) = qn(p,i,j,k) - dq1
|
||||
qb(p,it,js,ks) = qn(p,i,j,k) + dq2
|
||||
@ -6114,43 +6119,50 @@ module boundaries
|
||||
!
|
||||
do p = 1, nv
|
||||
|
||||
! calculate limited derivatives in all directions
|
||||
!
|
||||
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
|
||||
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
|
||||
dq(1) = limiter_prol(0.5d+00, dql, dqr)
|
||||
|
||||
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
|
||||
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
|
||||
dq(2) = limiter_prol(0.5d+00, dql, dqr)
|
||||
|
||||
dq(1) = 1.25d-01 * (qn(p,i+1,j,k) - qn(p,i-1,j,k))
|
||||
dq(2) = 1.25d-01 * (qn(p,i,j+1,k) - qn(p,i,j-1,k))
|
||||
#if NDIMS == 3
|
||||
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
|
||||
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
|
||||
dq(3) = limiter_prol(0.5d+00, dql, dqr)
|
||||
dq(3) = 1.25d-01 * (qn(p,i,j,k+1) - qn(p,i,j,k-1))
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (positive(p) .and. qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) then
|
||||
|
||||
if (qn(p,i,j,k) <= 0.0d+00) then
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
|
||||
dq(:) = 0.0d+00
|
||||
else
|
||||
|
||||
! calculate limited derivatives in all directions
|
||||
!
|
||||
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
|
||||
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
|
||||
dq(1) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
|
||||
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
|
||||
dq(2) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
#if NDIMS == 3
|
||||
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
|
||||
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
|
||||
dq(3) = limiter_prol(2.5d-01, dql, dqr)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (positive(p)) then
|
||||
if (qn(p,i,j,k) > 0.0d+00) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
|
||||
#if NDIMS == 2
|
||||
! calculate the derivative terms
|
||||
! calculate derivative terms
|
||||
!
|
||||
dq1 = dq(1) + dq(2)
|
||||
dq2 = dq(1) - dq(2)
|
||||
|
||||
! prolong the edge region to the output array
|
||||
! prolong the edge region
|
||||
!
|
||||
qb(p,is,js,k ) = qn(p,i,j,k) - dq1
|
||||
qb(p,it,js,k ) = qn(p,i,j,k) + dq2
|
||||
@ -6158,14 +6170,14 @@ module boundaries
|
||||
qb(p,it,jt,k ) = qn(p,i,j,k) + dq1
|
||||
#endif /* NDIMS == 2 */
|
||||
#if NDIMS == 3
|
||||
! calculate the derivative terms
|
||||
! calculate derivative terms
|
||||
!
|
||||
dq1 = dq(1) + dq(2) + dq(3)
|
||||
dq2 = dq(1) - dq(2) - dq(3)
|
||||
dq3 = dq(1) - dq(2) + dq(3)
|
||||
dq4 = dq(1) + dq(2) - dq(3)
|
||||
|
||||
! prolong the edge region to the output array
|
||||
! prolong the edge region
|
||||
!
|
||||
qb(p,is,js,ks) = qn(p,i,j,k) - dq1
|
||||
qb(p,it,js,ks) = qn(p,i,j,k) + dq2
|
||||
@ -6177,7 +6189,7 @@ module boundaries
|
||||
qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
end do ! q = 1, nv
|
||||
end do ! p = 1, nv
|
||||
|
||||
end do ! i = il, iu
|
||||
end do ! j = jl, ju
|
||||
@ -6391,43 +6403,50 @@ module boundaries
|
||||
!
|
||||
do p = 1, nv
|
||||
|
||||
! calculate limited derivatives in all directions
|
||||
!
|
||||
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
|
||||
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
|
||||
dq(1) = limiter_prol(0.5d+00, dql, dqr)
|
||||
|
||||
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
|
||||
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
|
||||
dq(2) = limiter_prol(0.5d+00, dql, dqr)
|
||||
|
||||
dq(1) = 1.25d-01 * (qn(p,i+1,j,k) - qn(p,i-1,j,k))
|
||||
dq(2) = 1.25d-01 * (qn(p,i,j+1,k) - qn(p,i,j-1,k))
|
||||
#if NDIMS == 3
|
||||
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
|
||||
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
|
||||
dq(3) = limiter_prol(0.5d+00, dql, dqr)
|
||||
dq(3) = 1.25d-01 * (qn(p,i,j,k+1) - qn(p,i,j,k-1))
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (positive(p) .and. qn(p,i,j,k) <= sum(abs(dq(1:NDIMS)))) then
|
||||
|
||||
if (qn(p,i,j,k) <= 0.0d+00) then
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
|
||||
dq(:) = 0.0d+00
|
||||
else
|
||||
|
||||
! calculate limited derivatives in all directions
|
||||
!
|
||||
dql = qn(p,i ,j,k) - qn(p,im1,j,k)
|
||||
dqr = qn(p,ip1,j,k) - qn(p,i ,j,k)
|
||||
dq(1) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
dql = qn(p,i,j ,k) - qn(p,i,jm1,k)
|
||||
dqr = qn(p,i,jp1,k) - qn(p,i,j ,k)
|
||||
dq(2) = limiter_prol(2.5d-01, dql, dqr)
|
||||
|
||||
#if NDIMS == 3
|
||||
dql = qn(p,i,j,k ) - qn(p,i,j,km1)
|
||||
dqr = qn(p,i,j,kp1) - qn(p,i,j,k )
|
||||
dq(3) = limiter_prol(2.5d-01, dql, dqr)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (positive(p)) then
|
||||
if (qn(p,i,j,k) > 0.0d+00) then
|
||||
do while (qn(p,i,j,k) <= sum(abs(dq(1:NDIMS))))
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
end do
|
||||
else
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
dq(:) = 0.0d+00
|
||||
end if
|
||||
end if
|
||||
|
||||
dq(:) = 0.5d+00 * dq(:)
|
||||
|
||||
#if NDIMS == 2
|
||||
! calculate the derivative terms
|
||||
! calculate derivative terms
|
||||
!
|
||||
dq1 = dq(1) + dq(2)
|
||||
dq2 = dq(1) - dq(2)
|
||||
|
||||
! prolong the corner region to the output array
|
||||
! prolong the corner region
|
||||
!
|
||||
qb(p,is,js,k ) = qn(p,i,j,k) - dq1
|
||||
qb(p,it,js,k ) = qn(p,i,j,k) + dq2
|
||||
@ -6435,14 +6454,14 @@ module boundaries
|
||||
qb(p,it,jt,k ) = qn(p,i,j,k) + dq1
|
||||
#endif /* NDIMS == 2 */
|
||||
#if NDIMS == 3
|
||||
! calculate the derivative terms
|
||||
! calculate derivative terms
|
||||
!
|
||||
dq1 = dq(1) + dq(2) + dq(3)
|
||||
dq2 = dq(1) - dq(2) - dq(3)
|
||||
dq3 = dq(1) - dq(2) + dq(3)
|
||||
dq4 = dq(1) + dq(2) - dq(3)
|
||||
|
||||
! prolong the corner region to the output array
|
||||
! prolong the corner region
|
||||
!
|
||||
qb(p,is,js,ks) = qn(p,i,j,k) - dq1
|
||||
qb(p,it,js,ks) = qn(p,i,j,k) + dq2
|
||||
@ -6454,7 +6473,7 @@ module boundaries
|
||||
qb(p,it,jt,kt) = qn(p,i,j,k) + dq1
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
end do ! q = 1, nv
|
||||
end do ! p = 1, nv
|
||||
|
||||
end do ! i = il, iu
|
||||
end do ! j = jl, ju
|
||||
|
||||
@ -45,7 +45,7 @@ program amun
|
||||
use domains , only : initialize_domains, finalize_domains
|
||||
use equations , only : initialize_equations, finalize_equations
|
||||
use equations , only : print_equations
|
||||
use equations , only : nv, nf
|
||||
use equations , only : nv, nf, errors
|
||||
use evolution , only : initialize_evolution, finalize_evolution
|
||||
use evolution , only : print_evolution
|
||||
use evolution , only : advance, new_time_step
|
||||
@ -118,6 +118,7 @@ program amun
|
||||
real(kind=8) :: zmin = 0.0d+00, zmax = 1.0d+00
|
||||
real(kind=8) :: tmax = 0.0d+00, trun = 9.999d+03, tsav = 3.0d+01
|
||||
real(kind=8) :: dtnext = 0.0d+00
|
||||
real(kind=8) :: maxerr = 0.0d+00
|
||||
|
||||
! timer indices
|
||||
!
|
||||
@ -613,16 +614,16 @@ program amun
|
||||
!
|
||||
write(*,*)
|
||||
write(*,"(1x,a)" ) "Evolving the system:"
|
||||
write(*,'(4x,a4,5x,a4,11x,a2,12x,a6,7x,a3)') 'step', 'time', 'dt' &
|
||||
, 'blocks', 'ETA'
|
||||
write(*,"(4x,'step',5x,'time',11x,'timestep',7x,'error',5x," // &
|
||||
"'blocks',7x,'ETA')")
|
||||
#ifdef __INTEL_COMPILER
|
||||
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1,$)') &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
|
||||
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1,$)") &
|
||||
step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#else /* __INTEL_COMPILER */
|
||||
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1)',advance="no") &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
|
||||
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1)",advance="no") &
|
||||
step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#endif /* __INTEL_COMPILER */
|
||||
|
||||
end if
|
||||
@ -652,6 +653,7 @@ program amun
|
||||
time = time + dt
|
||||
step = step + 1
|
||||
nsteps = nsteps + 1
|
||||
maxerr = maxval(errors)
|
||||
|
||||
! get current time in seconds
|
||||
!
|
||||
@ -706,13 +708,13 @@ program amun
|
||||
es = es - 60 * em
|
||||
|
||||
#ifdef __INTEL_COMPILER
|
||||
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1,$)') &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
|
||||
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1,$)") &
|
||||
step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#else /* __INTEL_COMPILER */
|
||||
write(*,'(i8,2(1x,1es14.6),2x,i8,2x,1i4.1,"d",1i2.2,"h",1i2.2,"m"' // &
|
||||
',1i2.2,"s",15x,a1)',advance="no") &
|
||||
step, time, dt, get_nleafs(), ed, eh, em, es, char(13)
|
||||
write(*,"(i8,2(1x,1es14.6),1x,1es10.2,2x,i8,2x," // &
|
||||
"1i4.1,'d',1i2.2,'h',1i2.2,'m',1i2.2,'s',15x,a1)",advance="no")&
|
||||
step, time, dt, maxerr, get_nleafs(), ed, eh, em, es, char(13)
|
||||
#endif /* __INTEL_COMPILER */
|
||||
|
||||
! update the timestamp
|
||||
@ -771,7 +773,7 @@ program amun
|
||||
"Problem finalizing FORCING module!"
|
||||
end if
|
||||
2100 continue
|
||||
call finalize_evolution(status)
|
||||
call finalize_evolution(master, status)
|
||||
if (check_status(status /= 0) .and. master) then
|
||||
write(error_unit,"('[AMUN::program]: ', a)") &
|
||||
"Problem finalizing EVOLUTION module!"
|
||||
|
||||
@ -221,6 +221,10 @@ module equations
|
||||
!
|
||||
logical, dimension(:), allocatable :: positive
|
||||
|
||||
! the array of maximum errors for conservative variables
|
||||
!
|
||||
real(kind=8), dimension(:), allocatable :: errors
|
||||
|
||||
! by default everything is private
|
||||
!
|
||||
private
|
||||
@ -250,7 +254,7 @@ module equations
|
||||
public :: eqsys, eos
|
||||
public :: pvars, cvars
|
||||
public :: qpbnd
|
||||
public :: ivars, positive
|
||||
public :: ivars, positive, errors
|
||||
|
||||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||||
!
|
||||
@ -1060,6 +1064,12 @@ module equations
|
||||
|
||||
end if
|
||||
|
||||
! allocate an array for errors
|
||||
!
|
||||
allocate(errors(nf), stat = status)
|
||||
|
||||
if (status == 0) errors(:) = 0.0d+00
|
||||
|
||||
! proceed if everything is fine
|
||||
!
|
||||
if (status == 0) then
|
||||
@ -1252,6 +1262,10 @@ module equations
|
||||
!
|
||||
if (allocated(positive)) deallocate(positive, stat = status)
|
||||
|
||||
! deallocate the array for errors
|
||||
!
|
||||
if (allocated(errors)) deallocate(errors, stat = status)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for module initialization/finalization
|
||||
!
|
||||
|
||||
@ -75,6 +75,21 @@ module evolution
|
||||
real(kind=8) , save :: time = 0.0d+00
|
||||
real(kind=8) , save :: dt = 1.0d+00
|
||||
real(kind=8) , save :: dtn = 1.0d+00
|
||||
real(kind=8) , save :: dte = 0.0d+00
|
||||
|
||||
! the absolute and relative tolerances, the maximum number of passes for
|
||||
! the adaptive step, the total number of integration iterations,
|
||||
! the number of rejected iterations
|
||||
!
|
||||
real(kind=8) , save :: atol = 1.0d-04
|
||||
real(kind=8) , save :: rtol = 1.0d-04
|
||||
integer , save :: mrej = 5
|
||||
integer , save :: niterations = 0
|
||||
integer , save :: nrejections = 0
|
||||
|
||||
! errors in three recent steps
|
||||
!
|
||||
real(kind=8), dimension(3), save :: errs
|
||||
|
||||
! the increment array
|
||||
!
|
||||
@ -91,7 +106,8 @@ module evolution
|
||||
|
||||
! declare public variables
|
||||
!
|
||||
public :: step, time, dt, dtn, cfl, glm_alpha
|
||||
public :: step, time, dt, dtn, dte, cfl, glm_alpha
|
||||
public :: atol, rtol, mrej, niterations, nrejections, errs
|
||||
|
||||
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||||
!
|
||||
@ -173,6 +189,13 @@ module evolution
|
||||
call get_parameter("cfl" , cfl )
|
||||
call get_parameter("glm_alpha" , glm_alpha )
|
||||
|
||||
! get the absolute and relative tolerances, the maximum number of passes for
|
||||
! the adaptive step
|
||||
!
|
||||
call get_parameter("absolute_tolerance", atol)
|
||||
call get_parameter("relative_tolerance", rtol)
|
||||
call get_parameter("maximum_rejections", mrej)
|
||||
|
||||
! select the integration method, check the correctness of the integration
|
||||
! parameters and adjust the CFL coefficient if necessary
|
||||
!
|
||||
@ -192,7 +215,7 @@ module evolution
|
||||
|
||||
stages = max(2, min(9, stages))
|
||||
write(name_int, "('2nd order SSPRK(',i0,',2)')") stages
|
||||
evolve => evolve_ssprk2
|
||||
evolve => evolve_ssprk2_m
|
||||
cfl = (stages - 1) * cfl
|
||||
|
||||
case ("rk3", "RK3")
|
||||
@ -270,6 +293,10 @@ module evolution
|
||||
|
||||
end if ! status
|
||||
|
||||
! reset recent error history
|
||||
!
|
||||
errs = 1.0d+00
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for module initialization/finalization
|
||||
!
|
||||
@ -293,10 +320,11 @@ module evolution
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine finalize_evolution(status)
|
||||
subroutine finalize_evolution(verbose, status)
|
||||
|
||||
! include external procedures
|
||||
!
|
||||
use helpers , only : print_section, print_parameter
|
||||
use iso_fortran_env, only : error_unit
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -305,6 +333,7 @@ module evolution
|
||||
|
||||
! subroutine arguments
|
||||
!
|
||||
logical, intent(in) :: verbose
|
||||
integer, intent(out) :: status
|
||||
|
||||
! local parameters
|
||||
@ -338,6 +367,14 @@ module evolution
|
||||
end if
|
||||
end if
|
||||
|
||||
! print integration summary
|
||||
!
|
||||
if (niterations > 0) then
|
||||
call print_section(verbose, "Integration summary")
|
||||
call print_parameter(verbose, "Number of iterations", niterations)
|
||||
call print_parameter(verbose, "Number of rejections", nrejections)
|
||||
end if
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for module initialization/finalization
|
||||
!
|
||||
@ -381,11 +418,14 @@ module evolution
|
||||
if (verbose) then
|
||||
|
||||
call print_section(verbose, "Evolution")
|
||||
call print_parameter(verbose, "time advance" , name_int)
|
||||
call print_parameter(verbose, "CFL coefficient" , cfl )
|
||||
call print_parameter(verbose, "time advance", name_int)
|
||||
call print_parameter(verbose, "CFL coefficient", cfl)
|
||||
if (magnetized) then
|
||||
call print_parameter(verbose, "GLM alpha coefficient", glm_alpha)
|
||||
end if
|
||||
call print_parameter(verbose, "absolute tolerance", atol)
|
||||
call print_parameter(verbose, "relative tolerance", rtol)
|
||||
call print_parameter(verbose, "maximum rejections", mrej)
|
||||
|
||||
end if
|
||||
|
||||
@ -609,7 +649,8 @@ module evolution
|
||||
|
||||
! round the time
|
||||
!
|
||||
if (dtnext > 0.0d+00) dt = min(dtn, dtnext)
|
||||
if (dte <= 0.0d+00) dte = dtn
|
||||
if (dtnext > 0.0d+00) dt = min(dtn, dte, dtnext)
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for new time step estimation
|
||||
@ -879,8 +920,8 @@ module evolution
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
! subroutine EVOLVE_SSPRK2:
|
||||
! ------------------------
|
||||
! subroutine EVOLVE_SSPRK2_M:
|
||||
! --------------------------
|
||||
!
|
||||
! Subroutine advances the solution by one time step using the 2nd order
|
||||
! m-stage Strong Stability Preserving Runge-Kutta time integration method.
|
||||
@ -896,180 +937,207 @@ module evolution
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
subroutine evolve_ssprk2()
|
||||
subroutine evolve_ssprk2_m()
|
||||
|
||||
! include external variables
|
||||
!
|
||||
use blocks , only : block_data, list_data
|
||||
use equations, only : ibp
|
||||
use sources , only : update_sources
|
||||
use blocks , only : block_data, list_data, get_dblocks
|
||||
use coordinates, only : nc => ncells, nb, ne
|
||||
use equations , only : errors, nf, ibp
|
||||
#ifdef MPI
|
||||
use mpitools , only : reduce_maximum
|
||||
#endif /* MPI */
|
||||
use sources , only : update_sources
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: n
|
||||
real(kind=8) :: ds
|
||||
real(kind=8) :: tm, dtm
|
||||
logical :: test
|
||||
integer :: n, l, nrej
|
||||
real(kind=8) :: tm, dtm, ds, umax, utol, emax
|
||||
|
||||
! local saved variables
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
real(kind=8), save :: ft, fl, fr
|
||||
real(kind=8), save :: ft, fl, fr, gt, gl
|
||||
|
||||
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! prepare things which don't change later
|
||||
!
|
||||
if (first) then
|
||||
|
||||
! calculate integration coefficients
|
||||
!
|
||||
ft = 1.0d+00 / (stages - 1)
|
||||
fl = 1.0d+00 / stages
|
||||
fr = 1.0d+00 - fl
|
||||
|
||||
! update first flag
|
||||
!
|
||||
gt = fl - ft
|
||||
gl = fl
|
||||
|
||||
first = .false.
|
||||
|
||||
end if
|
||||
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = ft * dt
|
||||
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 */
|
||||
|
||||
!= 1st step: U(0) = U(n)
|
||||
!
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
test = .true.
|
||||
nrej = 0
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
do while(test)
|
||||
|
||||
! copy conservative array u0 to u1
|
||||
!
|
||||
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
|
||||
lerr(:,:,:,:,:) = 0.0d+00
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pdata%u => pdata%u1
|
||||
ds = ft * dt
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
!= 2nd step: U(i) = [1 + dt/(m-1) L] U(i-1), for i = 1, ..., m-1
|
||||
!
|
||||
! integrate the intermediate steps
|
||||
!
|
||||
do n = 1, stages - 1
|
||||
|
||||
! prepare times
|
||||
!
|
||||
tm = time + n * ds
|
||||
dtm = ds
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
!= 1st step: U(1) = U(n), U(2) = U(n)
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_increment(pdata)
|
||||
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
|
||||
pdata%u2(:,:,:,:) = pdata%u0(:,:,:,:)
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
pdata%u => pdata%u1
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1
|
||||
!
|
||||
call update_variables(tm, dtm)
|
||||
do n = 1, stages - 1
|
||||
|
||||
end do ! n = 1, stages - 1
|
||||
tm = time + n * ds
|
||||
dtm = ds
|
||||
|
||||
!= the final step: U(n+1) = 1/m U(0) + (m-1)/m [1 + dt/(m-1) L] U(m-1)
|
||||
call update_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
call update_increment(pdata)
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! n = 1, stages - 1
|
||||
|
||||
!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)]
|
||||
!
|
||||
! prepare times
|
||||
tm = time + dt
|
||||
dtm = fl * dt
|
||||
|
||||
call update_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
call update_increment(pdata)
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
pdata%u1(:,:,:,:) = fr * pdata%u1(:,:,:,:) &
|
||||
+ fl * (pdata%u2(:,:,:,:) + dt * du(:,:,:,:))
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
! find umax
|
||||
!
|
||||
umax = 0.0d+00
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), &
|
||||
maxval(abs(pdata%u1(:,:,:,:))))
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
! get the maximum error of each variable
|
||||
!
|
||||
do l = 1, nf
|
||||
errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
|
||||
end do
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(umax)
|
||||
call reduce_maximum(errors)
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate tolerance and time step
|
||||
!
|
||||
utol = atol + rtol * umax
|
||||
emax = maxval(errors) / utol
|
||||
|
||||
if (emax <= 1.0d+00 .or. nrej >= mrej) then
|
||||
test = .false.
|
||||
|
||||
errs(3) = errs(2)
|
||||
errs(2) = errs(1)
|
||||
errs(1) = emax
|
||||
|
||||
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) &
|
||||
* errs(3)**(-0.025d+00)
|
||||
else
|
||||
errs(1) = emax
|
||||
|
||||
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) &
|
||||
* errs(3)**(-0.025d+00)
|
||||
dt = dte
|
||||
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
nrejections = nrejections + 1
|
||||
end if
|
||||
|
||||
niterations = niterations + 1
|
||||
|
||||
end do
|
||||
|
||||
if (allocated(lerr)) deallocate(lerr)
|
||||
|
||||
!= final step: U(n+1) = U(1)
|
||||
!
|
||||
tm = time + dt
|
||||
dtm = fl * dt
|
||||
dtm = ft * dt
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_increment(pdata)
|
||||
pdata%u0(:,:,:,:) = pdata%u1(:,:,:,:)
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(:,:,:,:) = fl * pdata%u0(:,:,:,:) &
|
||||
+ fr * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:))
|
||||
|
||||
! update the conservative variable pointer
|
||||
!
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
|
||||
|
||||
! assign pointer to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
#ifdef PROFILE
|
||||
@ -1080,7 +1148,7 @@ module evolution
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine evolve_ssprk2
|
||||
end subroutine evolve_ssprk2_m
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
@ -1852,42 +1920,34 @@ module evolution
|
||||
!
|
||||
subroutine evolve_ssprk3_m()
|
||||
|
||||
! include external variables
|
||||
!
|
||||
use blocks , only : block_data, list_data
|
||||
use equations, only : ibp
|
||||
use sources , only : update_sources
|
||||
use blocks , only : block_data, list_data, get_dblocks
|
||||
use coordinates, only : nc => ncells, nb, ne
|
||||
use equations , only : errors, nf, ibp
|
||||
#ifdef MPI
|
||||
use mpitools , only : reduce_maximum
|
||||
#endif /* MPI */
|
||||
use sources , only : update_sources
|
||||
|
||||
! local variables are not implicit by default
|
||||
!
|
||||
implicit none
|
||||
|
||||
! local pointers
|
||||
!
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
! local saved variables
|
||||
!
|
||||
logical , save :: first = .true.
|
||||
integer , save :: n, n1, n2, n3, n4
|
||||
real(kind=8), save :: r, f1, f2
|
||||
real(kind=8), save :: r, f1, f2, g1, g2
|
||||
|
||||
! local variables
|
||||
!
|
||||
integer :: i
|
||||
real(kind=8) :: ds
|
||||
real(kind=8) :: tm, dtm
|
||||
logical :: test
|
||||
integer :: i, l, nrej
|
||||
real(kind=8) :: tm, dtm, ds, umax, utol, emax
|
||||
|
||||
real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
#ifdef PROFILE
|
||||
! start accounting time for one step update
|
||||
!
|
||||
call start_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
! prepare coefficients
|
||||
!
|
||||
if (first) then
|
||||
|
||||
n = int(sqrt(1.0d+00 * stages))
|
||||
@ -1901,226 +1961,262 @@ module evolution
|
||||
f2 = (1.0d+00 * (n - 1)) / r
|
||||
r = 1.0d+00 * (stages - n)
|
||||
|
||||
! update first flag
|
||||
!
|
||||
g1 = 1.0d+00 / ( stages - n) - 1.0d+00 / stages
|
||||
g2 = 1.0d+00 / (2 * stages - n) - 1.0d+00 / stages
|
||||
|
||||
first = .false.
|
||||
|
||||
end if
|
||||
|
||||
!= 1st step: U(1) = [1 + dt/r) L] U(1), for i = 1, ..., n1
|
||||
!
|
||||
! calculate the fractional time step
|
||||
!
|
||||
ds = dt / r
|
||||
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 */
|
||||
|
||||
! integrate the intermediate steps
|
||||
!
|
||||
do i = 1, n1
|
||||
test = .true.
|
||||
nrej = 0
|
||||
|
||||
! prepare times
|
||||
!
|
||||
tm = time + i * ds
|
||||
dtm = ds
|
||||
do while(test)
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
lerr(:,:,:,:,:) = 0.0d+00
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
ds = dt / r
|
||||
|
||||
!= 1st step: U(1) = U(n)
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_increment(pdata)
|
||||
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
pdata%u => pdata%u1
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
! assign pdata to the next block
|
||||
!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2
|
||||
!
|
||||
do i = 1, n1
|
||||
|
||||
tm = time + i * ds
|
||||
dtm = ds
|
||||
|
||||
call update_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
call update_increment(pdata)
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! n = 1, n1
|
||||
|
||||
!= 3rd step: U(2) = U(1)
|
||||
!
|
||||
! iterate over all data blocks
|
||||
!
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
pdata%u2(:,:,:,:) = pdata%u1(:,:,:,:)
|
||||
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!= 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1
|
||||
!
|
||||
call update_variables(tm, dtm)
|
||||
do i = n2, n3
|
||||
|
||||
end do ! n = 1, n1
|
||||
tm = time + i * ds
|
||||
dtm = ds
|
||||
|
||||
!= 2nd step: U(2) = U(1)
|
||||
!
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
call update_fluxes()
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u1(:,:,:,:) = pdata%u0(:,:,:,:)
|
||||
call update_increment(pdata)
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
|
||||
|
||||
end do ! over data blocks
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
|
||||
!= 3rd step: U(1) = [1 + dt/r) L] U(1), for i = n2, ..., n3
|
||||
!
|
||||
! integrate the intermediate steps
|
||||
!
|
||||
do i = n2, n3
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
! prepare times
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! i = n2, n3
|
||||
|
||||
!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1)
|
||||
!
|
||||
tm = time + i * ds
|
||||
tm = time + dt
|
||||
dtm = ds
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_increment(pdata)
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
|
||||
pdata%u1(:,:,:,:) = f1 * pdata%u2(:,:,:,:) &
|
||||
+ f2 * (pdata%u1(:,:,:,:) + ds * du(:,:,:,:))
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! i = n2, n3
|
||||
|
||||
!= 4th step: U(1) = n * U(2) + (n - 1) * [1 + dt/r L] U(1)) / (2 * n - 1)
|
||||
!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m
|
||||
!
|
||||
! prepare times
|
||||
do i = n4, stages
|
||||
|
||||
tm = time + (i - n) * ds
|
||||
dtm = ds
|
||||
|
||||
call update_fluxes()
|
||||
|
||||
l = 1
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
|
||||
call update_increment(pdata)
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
pdata%u1(:,:,:,:) = pdata%u1(:,:,:,:) + ds * du(:,:,:,:)
|
||||
|
||||
#if NDIMS == 3
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * du(:,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
l = l + 1
|
||||
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
if (i < stages) call update_variables(tm, dtm)
|
||||
|
||||
end do ! i = n4, stages
|
||||
|
||||
! find umax
|
||||
!
|
||||
umax = 0.0d+00
|
||||
pdata => list_data
|
||||
do while (associated(pdata))
|
||||
umax = max(umax, maxval(abs(pdata%u0(:,:,:,:))), &
|
||||
maxval(abs(pdata%u1(:,:,:,:))))
|
||||
pdata => pdata%next
|
||||
end do ! over data blocks
|
||||
|
||||
! get the maximum error of each variable
|
||||
!
|
||||
do l = 1, nf
|
||||
errors(l) = dt * maxval(abs(lerr(:,l,:,:,:)))
|
||||
end do
|
||||
|
||||
#ifdef MPI
|
||||
call reduce_maximum(umax)
|
||||
call reduce_maximum(errors)
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate tolerance and time step
|
||||
!
|
||||
utol = atol + rtol * umax
|
||||
emax = maxval(errors) / utol
|
||||
|
||||
if (emax <= 1.0d+00 .or. nrej >= mrej) then
|
||||
test = .false.
|
||||
|
||||
errs(3) = errs(2)
|
||||
errs(2) = errs(1)
|
||||
errs(1) = emax
|
||||
|
||||
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) &
|
||||
* errs(3)**(-0.025d+00)
|
||||
else
|
||||
errs(1) = emax
|
||||
|
||||
dte = dt * errs(1)**(-0.2d+00) * errs(2)**(0.0505d+00) &
|
||||
* errs(3)**(-0.025d+00)
|
||||
dt = dte
|
||||
|
||||
nrej = nrej + 1 ! rejection count in the current step
|
||||
nrejections = nrejections + 1
|
||||
end if
|
||||
|
||||
niterations = niterations + 1
|
||||
|
||||
end do
|
||||
|
||||
if (allocated(lerr)) deallocate(lerr)
|
||||
|
||||
!= final step: U(n+1) = U(1)
|
||||
!
|
||||
tm = time + dt
|
||||
dtm = ds
|
||||
dtm = dt / r
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
!
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_increment(pdata)
|
||||
pdata%u0(:,:,:,:) = pdata%u1(:,:,:,:)
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
pdata%u => pdata%u0
|
||||
|
||||
! update the final solution
|
||||
!
|
||||
pdata%u0(:,:,:,:) = f1 * pdata%u1(:,:,:,:) + f2 * (pdata%u0(:,:,:,:) &
|
||||
+ ds * du(:,:,:,:))
|
||||
|
||||
! update ψ with its source term
|
||||
!
|
||||
if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
!= the final step: U(1) = [1 + dt/r) L] U(1), for i = n4, ..., m, U(n+1) = U(1)
|
||||
!
|
||||
! integrate the intermediate steps
|
||||
!
|
||||
do i = n4, stages
|
||||
|
||||
! prepare times
|
||||
!
|
||||
tm = time + (i - n) * ds
|
||||
dtm = ds
|
||||
|
||||
! update fluxes
|
||||
!
|
||||
call update_fluxes()
|
||||
|
||||
! assign pdata with the first block on the data block list
|
||||
! update ψ with its source term in the last step
|
||||
!
|
||||
if (ibp > 0) then
|
||||
pdata => list_data
|
||||
|
||||
! iterate over all data blocks
|
||||
!
|
||||
do while (associated(pdata))
|
||||
|
||||
! calculate the variable increment
|
||||
!
|
||||
call update_increment(pdata)
|
||||
|
||||
! add the source terms
|
||||
!
|
||||
call update_sources(pdata, tm, dtm, du(:,:,:,:))
|
||||
|
||||
! update the intermediate solution
|
||||
!
|
||||
pdata%u0(:,:,:,:) = pdata%u0(:,:,:,:) + ds * du(:,:,:,:)
|
||||
|
||||
! assign pdata to the next block
|
||||
!
|
||||
pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:)
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! over data blocks
|
||||
end if
|
||||
|
||||
! update primitive variables
|
||||
!
|
||||
call update_variables(tm, dtm)
|
||||
|
||||
end do ! i = n4, stages
|
||||
|
||||
#ifdef PROFILE
|
||||
! stop accounting time for one step update
|
||||
!
|
||||
call stop_timer(imu)
|
||||
#endif /* PROFILE */
|
||||
|
||||
|
||||
@ -51,17 +51,23 @@ module integrals
|
||||
!
|
||||
! funit - a file handler to the integrals file;
|
||||
! sunit - a file handler to the statistics file;
|
||||
! eunit - a file handler to the errors file;
|
||||
! iintd - the number of steps between subsequent intervals storing;
|
||||
!
|
||||
integer(kind=4), save :: funit = 11
|
||||
integer(kind=4), save :: sunit = 12
|
||||
integer(kind=4), save :: runit = 13
|
||||
integer(kind=4), save :: eunit = 13
|
||||
integer(kind=4), save :: runit = 14
|
||||
integer(kind=4), save :: iintd = 1
|
||||
|
||||
! flag indicating if the files are actually written to disk
|
||||
!
|
||||
logical, save :: stored = .false.
|
||||
|
||||
! the format string
|
||||
!
|
||||
character(len=32), save :: efmt
|
||||
|
||||
! by default everything is private
|
||||
!
|
||||
private
|
||||
@ -100,6 +106,7 @@ module integrals
|
||||
|
||||
! import external variables and subroutines
|
||||
!
|
||||
use equations , only : pvars, nf
|
||||
use parameters, only : get_parameter
|
||||
|
||||
! local variables are not implicit by default
|
||||
@ -114,7 +121,7 @@ module integrals
|
||||
|
||||
! local variables
|
||||
!
|
||||
character(len=32) :: fname, append
|
||||
character(len=32) :: fname, append, stmp
|
||||
logical :: flag
|
||||
!
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -240,6 +247,53 @@ module integrals
|
||||
, 'mean(Malf)', 'min(Malf)', 'max(Malf)'
|
||||
write(sunit,"('#')")
|
||||
|
||||
! depending on the append parameter, choose the right file initialization for
|
||||
! the file to store integration errors
|
||||
!
|
||||
append = "off"
|
||||
call get_parameter("errors_append", append)
|
||||
select case(trim(append))
|
||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||
write(fname, "('errors.dat')")
|
||||
inquire(file = fname, exist = flag)
|
||||
case default
|
||||
write(fname, "('errors_',i2.2,'.dat')") irun
|
||||
flag = .false.
|
||||
end select
|
||||
|
||||
! check if the file exists; if not, create a new one, otherwise move to the end
|
||||
!
|
||||
if (flag .and. irun > 1) then
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit = eunit, file = fname, form = 'formatted', status = 'old' &
|
||||
, position = 'append', buffered = 'yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit = eunit, file = fname, form = 'formatted', status = 'old' &
|
||||
, position = 'append')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
write(funit,"('#')")
|
||||
else
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit = eunit, file = fname, form = 'formatted' &
|
||||
, status = 'replace', buffered = 'yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit = eunit, file = fname, form = 'formatted' &
|
||||
, status = 'replace')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
end if
|
||||
|
||||
! write the integral file header
|
||||
!
|
||||
write(stmp,"(i9)") nf + 4
|
||||
efmt = "('#',a8," // trim(adjustl(stmp)) // "(1x,a18))"
|
||||
write(eunit,efmt) 'step', 'time', 'dt_cfl', 'dt_err', &
|
||||
'max_err', pvars(1:nf)
|
||||
write(eunit,"('#')")
|
||||
|
||||
! prepare format for errors
|
||||
!
|
||||
efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))"
|
||||
|
||||
! depending on the append parameter, choose the right file initialization for
|
||||
! the reconnection file
|
||||
!
|
||||
@ -336,6 +390,7 @@ module integrals
|
||||
if (stored) then
|
||||
close(funit)
|
||||
close(sunit)
|
||||
close(eunit)
|
||||
close(runit)
|
||||
end if
|
||||
|
||||
@ -378,7 +433,8 @@ module integrals
|
||||
use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
|
||||
use equations , only : ien, imx, imy, imz
|
||||
use equations , only : magnetized, adiabatic_index, csnd
|
||||
use evolution , only : step, time, dtn
|
||||
use equations , only : errors
|
||||
use evolution , only : step, time, dtn, dte
|
||||
use forcing , only : einj, rinj, arms
|
||||
#ifdef MPI
|
||||
use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum
|
||||
@ -903,6 +959,7 @@ module integrals
|
||||
, avarr(5), mnarr(5), mxarr(5) &
|
||||
, avarr(6), mnarr(6), mxarr(6) &
|
||||
, avarr(7), mnarr(7), mxarr(7)
|
||||
write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:)
|
||||
write(runit,"(i9, 11es18.8e3)") step, time, inarr(11:20)
|
||||
end if
|
||||
|
||||
|
||||
@ -1261,6 +1261,7 @@ module io
|
||||
use coordinates , only : nn => bcells, ncells
|
||||
use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax
|
||||
use evolution , only : step, time, dt, dtn
|
||||
use evolution , only : niterations, nrejections, errs
|
||||
use forcing , only : nmodes, fcoefs, einj
|
||||
use hash , only : xxh64
|
||||
use iso_fortran_env, only : error_unit
|
||||
@ -1412,6 +1413,16 @@ module io
|
||||
read(svalue, fmt = *) dt
|
||||
case('dtn')
|
||||
read(svalue, fmt = *) dtn
|
||||
case('niterations')
|
||||
read(svalue, fmt = *) niterations
|
||||
case('nrejections')
|
||||
read(svalue, fmt = *) nrejections
|
||||
case('errs(1)')
|
||||
read(svalue, fmt = *) errs(1)
|
||||
case('errs(2)')
|
||||
read(svalue, fmt = *) errs(2)
|
||||
case('errs(3)')
|
||||
read(svalue, fmt = *) errs(3)
|
||||
case('fields')
|
||||
il = index(line, 'digest="') + 8
|
||||
iu = index(line(il:), '"') + il - 2
|
||||
@ -2065,7 +2076,8 @@ module io
|
||||
#endif /* NDIMS == 3 */
|
||||
use coordinates , only : bdims => domain_base_dims
|
||||
use equations , only : eqsys, eos, nv
|
||||
use evolution , only : step, time, dt, dtn, cfl, glm_alpha
|
||||
use evolution , only : step, time, dt, dtn, cfl, glm_alpha, errs
|
||||
use evolution , only : atol, rtol, mrej, niterations, nrejections
|
||||
use forcing , only : nmodes, fcoefs, einj
|
||||
use iso_fortran_env, only : error_unit
|
||||
use mpitools , only : nprocs, nproc
|
||||
@ -2216,6 +2228,14 @@ module io
|
||||
call write_attribute_xml(lun, "dtn" , dtn)
|
||||
call write_attribute_xml(lun, "cfl" , cfl)
|
||||
call write_attribute_xml(lun, "glm_alpha", glm_alpha)
|
||||
call write_attribute_xml(lun, "absolute_tolerance", atol)
|
||||
call write_attribute_xml(lun, "relative_tolerance", rtol)
|
||||
call write_attribute_xml(lun, "maximum_rejections", mrej)
|
||||
call write_attribute_xml(lun, "niterations", niterations)
|
||||
call write_attribute_xml(lun, "nrejections", nrejections)
|
||||
call write_attribute_xml(lun, "errs(1)", errs(1))
|
||||
call write_attribute_xml(lun, "errs(2)", errs(2))
|
||||
call write_attribute_xml(lun, "errs(3)", errs(3))
|
||||
write(lun,"(a)") '</Evolution>'
|
||||
write(lun,"(a)") '<Forcing>'
|
||||
call write_attribute_xml(lun, "nmodes" , nmodes)
|
||||
@ -3990,7 +4010,8 @@ module io
|
||||
use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax
|
||||
use coordinates , only : periodic
|
||||
use equations , only : eqsys, eos, adiabatic_index, csnd
|
||||
use evolution , only : step, time, dt, dtn, cfl, glm_alpha
|
||||
use evolution , only : step, time, dt, dtn, cfl, glm_alpha, errs
|
||||
use evolution , only : atol, rtol, mrej, niterations, nrejections
|
||||
use forcing , only : nmodes, einj, fcoefs
|
||||
use hdf5 , only : hid_t
|
||||
use hdf5 , only : h5gcreate_f, h5gclose_f
|
||||
@ -4078,6 +4099,9 @@ module io
|
||||
call write_attribute(gid, 'step' , step )
|
||||
call write_attribute(gid, 'isnap' , isnap )
|
||||
call write_attribute(gid, 'periodic', per(:) )
|
||||
call write_attribute(gid, 'niterations', niterations)
|
||||
call write_attribute(gid, 'nrejections', nrejections)
|
||||
call write_attribute(gid, 'maximum_rejections', mrej)
|
||||
|
||||
! store the real attributes
|
||||
!
|
||||
@ -4098,6 +4122,11 @@ module io
|
||||
if (eos == 'iso') then
|
||||
call write_attribute(gid, 'sound_speed', csnd)
|
||||
end if
|
||||
call write_attribute(gid, 'absolute_tolerance', atol)
|
||||
call write_attribute(gid, 'relative_tolerance', rtol)
|
||||
call write_attribute(gid, 'errs(1)', errs(1))
|
||||
call write_attribute(gid, 'errs(2)', errs(2))
|
||||
call write_attribute(gid, 'errs(3)', errs(3))
|
||||
|
||||
! store the vector attributes
|
||||
!
|
||||
@ -4182,6 +4211,7 @@ module io
|
||||
use coordinates , only : ncells
|
||||
use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax
|
||||
use evolution , only : step, time, dt, dtn
|
||||
use evolution , only : niterations, nrejections, errs
|
||||
use forcing , only : nmodes, fcoefs
|
||||
use hdf5 , only : hid_t
|
||||
use hdf5 , only : h5gopen_f, h5gclose_f
|
||||
@ -4242,6 +4272,8 @@ module io
|
||||
call read_attribute(gid, 'nseeds' , lnseeds )
|
||||
call read_attribute(gid, 'step' , step )
|
||||
call read_attribute(gid, 'isnap' , isnap )
|
||||
call read_attribute(gid, 'niterations', niterations)
|
||||
call read_attribute(gid, 'nrejections', nrejections)
|
||||
|
||||
! restore double precision attributes
|
||||
!
|
||||
@ -4254,6 +4286,9 @@ module io
|
||||
call read_attribute(gid, 'time', time)
|
||||
call read_attribute(gid, 'dt' , dt )
|
||||
call read_attribute(gid, 'dtn' , dtn )
|
||||
call read_attribute(gid, 'errs(1)', errs(1))
|
||||
call read_attribute(gid, 'errs(2)', errs(2))
|
||||
call read_attribute(gid, 'errs(3)', errs(3))
|
||||
|
||||
! check the number of dimensions
|
||||
!
|
||||
|
||||
@ -1190,36 +1190,39 @@ module mesh
|
||||
|
||||
do p = 1, nv
|
||||
|
||||
dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k)
|
||||
dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k)
|
||||
du(1) = limiter_prol(0.5d+00, dul, dur)
|
||||
|
||||
dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k)
|
||||
dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k)
|
||||
du(2) = limiter_prol(0.5d+00, dul, dur)
|
||||
|
||||
du(1) = 1.25d-01 * (pdata%u(p,i+1,j,k) - pdata%u(p,i-1,j,k))
|
||||
du(2) = 1.25d-01 * (pdata%u(p,i,j+1,k) - pdata%u(p,i,j-1,k))
|
||||
#if NDIMS == 3
|
||||
dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1)
|
||||
dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k )
|
||||
du(3) = limiter_prol(0.5d+00, dul, dur)
|
||||
du(3) = 1.25d-01 * (pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k-1))
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (positive(p)) then
|
||||
if (positive(p) .and. pdata%u(p,i,j,k) < sum(abs(du(1:NDIMS)))) then
|
||||
if (pdata%u(p,i,j,k) > 0.0d+00) then
|
||||
|
||||
dul = pdata%u(p,i ,j,k) - pdata%u(p,i-1,j,k)
|
||||
dur = pdata%u(p,i+1,j,k) - pdata%u(p,i ,j,k)
|
||||
du(1) = limiter_prol(2.5d-01, dul, dur)
|
||||
|
||||
dul = pdata%u(p,i,j ,k) - pdata%u(p,i,j-1,k)
|
||||
dur = pdata%u(p,i,j+1,k) - pdata%u(p,i,j ,k)
|
||||
du(2) = limiter_prol(2.5d-01, dul, dur)
|
||||
|
||||
#if NDIMS == 3
|
||||
dul = pdata%u(p,i,j,k ) - pdata%u(p,i,j,k-1)
|
||||
dur = pdata%u(p,i,j,k+1) - pdata%u(p,i,j,k )
|
||||
du(3) = limiter_prol(2.5d-01, dul, dur)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
do while (pdata%u(p,i,j,k) <= sum(abs(du(1:NDIMS))))
|
||||
du(:) = 0.5d+00 * du(:)
|
||||
end do
|
||||
else
|
||||
write(error_unit,"('[',a,']: ',a,i9,a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive for block ID:" &
|
||||
, pmeta%id, " at ( ", i, j, k, " )"
|
||||
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc) &
|
||||
, "Positive variable is not positive at (", i, j, k, " )"
|
||||
du(:) = 0.0d+00
|
||||
status = 1
|
||||
end if
|
||||
end if
|
||||
|
||||
du(:) = 0.5d+00 * du(:)
|
||||
|
||||
#if NDIMS == 2
|
||||
du1 = du(1) + du(2)
|
||||
du2 = du(1) - du(2)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user