! $Id$
!
module Timestep
!
use Cparam
use Cdata
!
implicit none
!
private
!
include 'timestep.h'
!
! Parameters for adaptive time stepping
real, parameter :: safety = 0.9
real, parameter :: dt_decrease = -0.25
real, parameter :: dt_increase = -0.20
real :: errcon
!
contains
!
!***********************************************************************
subroutine initialize_timestep
endsubroutine initialize_timestep
!***********************************************************************
subroutine time_step(f,df,p)
!
! Runge-Kutta-Fehlberg accurate to 5th order
!
! 22-jun-06/tony: coded
!
use Cdata
use Messages
!! use Particles_main
!! use Shear, only: advance_shear
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
real :: errmax, tnew
real :: dt_temp, dt_next, dt_did
integer :: j,i
!
ldt=.false.
!
! General error condition
errcon = (5.0/safety)**(1.0/dt_increase)
!
if (itorder/=5) &
call fatal_error('time_step','itorder must be 5 for Runge-Kutta-Fehlberg')
!
!
! dt_beta_ts may be needed in other modules (like Dustdensity) for fixed dt
!
! if (.not. ldt) dt_beta_ts=dt*beta_ts
!
!
if (linterstellar.or.lshear.or.lparticles) &
call fatal_error("time_step", &
"Shear, interstallar and particles are not" // &
" yet supported by the adaptive rkf scheme")
!
lfirst=.true.
do i=1,10
! Do a Runge-Kutta step
call rkck(f, df, p, errmax)
! Step succeeded so exit
if (errmax <= safety) exit
! Step didn't succeed so decrease the time step
! print*,"Decreasing"
dt_temp = safety*dt*(errmax**dt_decrease)
! Don't decrease the time step by more than a factor of ten
dt = sign(max(abs(dt_temp), 0.1*abs(dt)), dt)
! New time
tnew = t+dt
if (tnew == t) then
! Guard against infinitesimal time steps
print*, 'WARNING: Timestep underflow in rkqs()'
endif
enddo
!
! print*,"errmax, errcon", errmax,errcon
if (errmax > errcon) then
! Increase the time step
dt_next = safety*dt*(errmax**dt_increase)
else
! But not by more than a factor of 5
dt_next = 5.0*dt
endif
!
! Time step that was actually performed
dt_did = dt
!
if (ip<=6) print*,'TIMESTEP: iproc,dt=',iproc_world,dt !(all have same dt?)
! Increase time
t = t+dt
! Time step to try next time
dt = dt_next
!
! Time evolution of grid variables
! (do this loop in pencils, for cache efficiency)
!
do j=1,mvar; do n=n1,n2; do m=m1,m2
f(l1:l2,m,n,j)=f(l1:l2,m,n,j)+df(l1:l2,m,n,j)
enddo; enddo; enddo
!
endsubroutine time_step
!***********************************************************************
subroutine rkck(f, df, p, errmax)
! Explicit fifth order Runge--Kutta--Fehlberg time stepping
use Cdata
use Mpicomm, only: mpiallreduce_max,MPI_COMM_WORLD
use Messages
use Equ
! RK parameters by Cash and Karp
real, parameter :: b21 = 0.2
real, parameter :: b31 = 0.075
real, parameter :: b32 = 0.225
real, parameter :: b41 = 0.3
real, parameter :: b42 = -0.9
real, parameter :: b43 = 1.2
real, parameter :: b51 = -11.0 / 54.0
real, parameter :: b52 = 2.5
real, parameter :: b53 = -70.0 / 27.0
real, parameter :: b54 = 35.0 / 27.0
real, parameter :: b61 = 1631.0 / 55296.0
real, parameter :: b62 = 175.0 / 512.0
real, parameter :: b63 = 575.0 / 13824.0
real, parameter :: b64 = 44275.0 / 110592.0
real, parameter :: b65 = 253.0 / 4096.0
real, parameter :: c1 = 37.0 / 378.0
real, parameter :: c2 = 0.0
real, parameter :: c3 = 250.0 / 621.0
real, parameter :: c4 = 125.0 / 594.0
real, parameter :: c5 = 0.0
real, parameter :: c6 = 512.0 / 1771.0
real, parameter :: dc1 = c1 - 2825.0 / 27648.0
real, parameter :: dc2 = c2 - 0.0
real, parameter :: dc3 = c3 - 18575.0 / 48384.0
real, parameter :: dc4 = c4 - 13525.0 / 55296.0
real, parameter :: dc5 = c5 - 277.0 / 14336.0
real, parameter :: dc6 = c6 - 0.25
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
real, dimension (mx,my,mz,mvar), intent(out) :: df
type (pencil_case), intent(inout) :: p
real, dimension(mx,my,mz,mfarray,5) :: k
real, save, dimension(mx,my,mz,mfarray) :: tmp1
real, dimension(mx,my,mz,mfarray) :: tmp2
real, dimension(nx) :: scal, err
real, intent(inout) :: errmax
real :: errmaxs
integer :: j
logical, save :: first_call=.true.
!
if (ny /= 1 .or. nz /= 1) then
call fatal_error("rkck", "timestep_rkf_1d only works for the 1D case")
endif
!
errmax=0.
!
if (first_call) then
! Initialize tmp1 to arbitrary value /= 0, so chemistry.f90 (which
! operates in the ghost zones where it most probably shouldn't)
! doesn't divide by 0.
tmp1 = real(0.577215664901532860606512090082402431042159335)
first_call = .false.
endif
!
tmp2(:,m1:m2,n1:n2,:) = 0.
call pde(f,tmp2,p)
k(:,m1:m2,n1:n2,:,1) = tmp2(:,m1:m2,n1:n2,:)
do j=1,mvar; do n=n1,n2; do m=m1,m2
k(l1:l2,m,n,j,1) = dt*k(l1:l2,m,n,j,1)
!
enddo; enddo; enddo
!
lfirst=.false.
!
tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
+ b21*k(:,m1:m2,n1:n2,:,1)
tmp2(:,m1:m2,n1:n2,:) = 0.
call pde(tmp1,tmp2,p)
k(:,m1:m2,n1:n2,:,2) = tmp2(:,m1:m2,n1:n2,:)
do j=1,mvar; do n=n1,n2; do m=m1,m2
k(l1:l2,m,n,j,2) = dt*k(l1:l2,m,n,j,2)
enddo; enddo; enddo
!
tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
+ b31*k(:,m1:m2,n1:n2,:,1) &
+ b32*k(:,m1:m2,n1:n2,:,2)
tmp2(:,m1:m2,n1:n2,:) = 0.
call pde(tmp1,tmp2,p)
k(:,m1:m2,n1:n2,:,3) = tmp2(:,m1:m2,n1:n2,:)
do j=1,mvar; do n=n1,n2; do m=m1,m2
k(l1:l2,m,n,j,3) = dt*k(l1:l2,m,n,j,3)
enddo; enddo; enddo
!
tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
+ b41*k(:,m1:m2,n1:n2,:,1) &
+ b42*k(:,m1:m2,n1:n2,:,2) &
+ b43*k(:,m1:m2,n1:n2,:,3)
tmp2(:,m1:m2,n1:n2,:) = 0.
call pde(tmp1,tmp2,p)
k(:,m1:m2,n1:n2,:,4) = tmp2(:,m1:m2,n1:n2,:)
do j=1,mvar; do n=n1,n2; do m=m1,m2
k(l1:l2,m,n,j,4) = dt*k(l1:l2,m,n,j,4)
enddo; enddo; enddo
!
tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
+ b51*k(:,m1:m2,n1:n2,:,1) &
+ b52*k(:,m1:m2,n1:n2,:,2) &
+ b53*k(:,m1:m2,n1:n2,:,3) &
+ b54*k(:,m1:m2,n1:n2,:,4)
tmp2(:,m1:m2,n1:n2,:) = 0.
call pde(tmp1,tmp2,p)
k(:,m1:m2,n1:n2,:,5) = tmp2(:,m1:m2,n1:n2,:)
do j=1,mvar; do n=n1,n2; do m=m1,m2
k(l1:l2,m,n,j,5) = dt*k(l1:l2,m,n,j,5)
enddo; enddo; enddo
!
errmaxs=0.
tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
+ b61*k(:,m1:m2,n1:n2,:,1) &
+ b62*k(:,m1:m2,n1:n2,:,2) &
+ b63*k(:,m1:m2,n1:n2,:,3) &
+ b64*k(:,m1:m2,n1:n2,:,4) &
+ b65*k(:,m1:m2,n1:n2,:,5)
df(:,m1:m2,n1:n2,:)=0.
call pde(tmp1,df,p)
do j=1,mvar; do n=n1,n2; do m=m1,m2
df(l1:l2,m,n,j) = dt*df(l1:l2,m,n,j)
!
err = dc1*k(l1:l2,m,n,j,1) + dc2*k(l1:l2,m,n,j,2) + &
dc3*k(l1:l2,m,n,j,3) + dc4*k(l1:l2,m,n,j,4) + &
dc5*k(l1:l2,m,n,j,5) + dc6*df(l1:l2,m,n,j)
!
df(l1:l2,m,n,j) = c1*k(l1:l2,m,n,j,1) + c2*k(l1:l2,m,n,j,2) + &
c3*k(l1:l2,m,n,j,3) + c4*k(l1:l2,m,n,j,4) + &
c5*k(l1:l2,m,n,j,5) + c6*df(l1:l2,m,n,j)
!
! Get the maximum error over the whole field
!
select case (timestep_scaling(j))
case ('per_var_err')
!
! Per variable error
!
scal= ( &
sqrt(f(l1:l2,m,n,1)**2+f(l1:l2,m,n,2)**2) + &
sqrt(k(l1:l2,m,n,1,1)**2 + k(l1:l2,m,n,2,1)**2) + &
1e-30)
errmaxs = max(maxval(abs(err/scal)),errmaxs)
!scal= ( &
! abs(f(l1:l2,m,n,j)) + abs(k(l1:l2,m,n,j,1)) + 1e-30)
!errmaxs = max(maxval(abs(err/scal)),errmaxs)
case ('cons_frac_err')
!
! Constant fractional error
!
errmaxs = max(maxval(abs(err/f(l1:l2,m,n,j))),errmaxs)
case ('cons_err')
!
! Constant error
!
scal = max(abs(f(l1:l2,n,m,j)), 1e-8)
errmaxs = max(maxval(abs(err/scal)),errmaxs)
!
case ('none')
!
! No error check
!
errmaxs = 0
!
endselect
!
enddo; enddo; enddo
!
! Divide your maximum error by the required accuracy
!
errmaxs=errmaxs/eps_rkf
!
call mpiallreduce_max(errmaxs,errmax,MPI_COMM_WORLD)
!
endsubroutine rkck
!***********************************************************************
subroutine pushpars2c(p_par)
use Messages, only: fatal_error
integer, parameter :: n_pars=0
integer(KIND=ikind8), dimension(:) :: p_par
call fatal_error('timestep_rkf_1d','alpha_ts, beta_ts not defined')
endsubroutine pushpars2c
!***********************************************************************
endmodule Timestep