! $Id$
!
! MODULE_DOC: This is a highly specified timestep module currently only working
! MODULE_DOC: together with the special module coronae.f90.
!
module Timestep
!
use Cparam
use Cdata
!
implicit none
!
private
!
include 'timestep.h'
!
contains
!***********************************************************************
subroutine initialize_timestep
!
! Coefficients for order 3.
!
use Messages, only: fatal_error
use General, only: itoa
!
if (itorder==3) then
!
! Coefficients for order 3. (use coefficients of Williamson (1980))
!
alpha_ts=(/ 0.0, -5/9.0 , -153/128.0, 0., 0. /)
beta_ts =(/ 1/3.0, 15/16.0, 8/15.0, 0., 0. /)
else
call fatal_error('initialize_timestep','Not implemented: itorder= '// &
trim(itoa(itorder)))
endif
endsubroutine initialize_timestep
!***********************************************************************
subroutine time_step(f,df,p)
!
! 2-apr-01/axel: coded
! 14-sep-01/axel: moved itorder to cdata
!
use Boundcond
use BorderProfiles, only: border_quenching
use Equ, only: pde
use Mpicomm, only: mpiallreduce_max
use Particles_main, only: particles_timestep_first, &
particles_timestep_second
use Shear, only: advance_shear
use Sub, only: shift_dt
use Energy
use Special
use Boundcond
use General, only: notanumber
!
real, dimension (mx,my,mz,mfarray) :: f,fsub,fm1
real, dimension (mx,my,mz) :: fm2, fj
real, dimension (mx,my,mz,mvar) :: df,dfsub,dfm1
type (pencil_case) :: p
real :: ds
real :: dt1, dt1_local, dt1_last=0.0
integer :: ienergy, Nsub, itRKL, s
real, dimension (nx) :: dt1_hcond_max
real, dimension (itorder) :: coeff_fm1, coeff_fm2, coeff_fsub, coeff_dfm1, coeff_dfsub
real :: dt1_energy_local,dt1_energy,dt_energy,dt_RKL,dt_sub_RKL
!
ienergy = ilnTT
!
! dt_beta_ts may be needed in other modules (like Dustdensity) for fixed dt.
!
if (.not. ldt) dt_beta_ts=dt*beta_ts
!
! Set up df and ds for each time sub.
!
do itsub=1,3
lfirst=(itsub==1)
llast=(itsub==3)
if (lfirst) then
df=0.0
ds=0.0
else
df=alpha_ts(itsub)*df !(could be subsumed into pde, but is dangerous!)
ds=alpha_ts(itsub)*ds
endif
!
! Change df according to the chosen physics modules.
!
call pde(f,df,p)
!
ds=ds+1.0
!
! If we are in the first time substep we need to calculate timestep dt.
! Only do it on the root processor, then broadcast dt to all others.
!
if (lfirst.and.ldt) then
dt1_local=maxval(dt1_max(1:nx))
! Timestep growth limiter
if (real(ddt) > 0.) dt1_local=max(dt1_local,dt1_last)
! Get global time step limite
call mpiallreduce_max(dt1_local,dt1,MPI_COMM_WORLD)
dt=1.0/dt1
! in pde(f,df,p) ghost cells of f-array are set
fsub = f
if (loutput_varn_at_exact_tsnap) call shift_dt(dt)
! Timestep growth limiter
if (ddt/=0.) dt1_last=dt1_local/ddt
endif
!
! Calculate dt_beta_ts.
!
if (ldt) dt_beta_ts=dt*beta_ts
if (ip<=6) print*, 'time_step: iproc, dt=', iproc_world, dt !(all have same dt?)
!
! Apply border quenching.
!
if (lborder_profiles) call border_quenching(f,df,dt_beta_ts(itsub))
!
! Time evolution of grid variables.
!
f(l1:l2,m1:m2,n1:n2,1:mvar) = f(l1:l2,m1:m2,n1:n2,1:mvar) + dt_beta_ts(itsub)*df(l1:l2,m1:m2,n1:n2,1:mvar)
!
if (lspecial) call special_after_timestep(f,df,dt_beta_ts(itsub)*ds,llast)
!
! Increase time.
!
t = t + dt_beta_ts(itsub)*ds
!
enddo
!--------------------------------------------------------------------------------*
! Now do the substeps for the energy (thermal conduction) only -----------------*
!--------------------------------------------------------------------------------*
!
if (sub_step_hcond) then
!
! initialized fsub
fsub(l1:l2,m1:m2,n1:n2,ienergy)=f(l1:l2,m1:m2,n1:n2,ienergy)
! update boundary is necessary for calculate time steps
call boundconds_x(fsub,ilnTT,ilnTT)
call initiate_isendrcv_bdry(fsub,ilnTT,ilnTT)
call finalize_isendrcv_bdry(fsub,ilnTT,ilnTT)
call boundconds_y(fsub,ilnTT,ilnTT)
call boundconds_z(fsub,ilnTT,ilnTT)
!
! This should be a better position to determin dt_energy, changes in lnTT
! by other terms are taken into account.
! Get time step for heat conduction in sub step
call calc_hcond_timestep(fsub,p,dt1_hcond_max)
dt1_energy_local=maxval(dt1_hcond_max(1:nx))
call mpiallreduce_max(dt1_energy_local,dt1_energy,MPI_COMM_WORLD)
!
dt_energy = 1d0/dt1_energy
! Set time step to the super-timestep
itRKL = itorder
dt_RKL = (itRKL*itRKL+itRKL-2d0)/4d0*dt_energy
!
! calc the number of sub iteration
Nsub = int(dt / dt_RKL)+1
dt_sub_RKL = dt / dble(Nsub)
!
! calc timestep coefficient
call RKL_coeff(itRKL, coeff_fm1, coeff_fm2, coeff_fsub, coeff_dfm1, coeff_dfsub)
!
! initalize fsub, fm1, fm2 for sub cycle with RKL steps
fm1=fsub
!
do j=1,Nsub
!
! one RKL step contains itRKL sub steps
!
do s=1,itRKL
if (s == 1) then
dfsub=0.
call pde_energy_only(fsub,dfsub,p,dt_sub_RKL)
fm1(l1:l2,m1:m2,n1:n2,ienergy)=fsub(l1:l2,m1:m2,n1:n2,ienergy) &
+coeff_dfm1(s)*dt_sub_RKL*dfsub(l1:l2,m1:m2,n1:n2,ienergy)
fm2(l1:l2,m1:m2,n1:n2) = fsub(l1:l2,m1:m2,n1:n2,ienergy)
else
dfm1=0.
call pde_energy_only(fm1,dfm1,p,dt_sub_RKL)
fj(l1:l2,m1:m2,n1:n2)=coeff_fm1(s) * fm1(l1:l2,m1:m2,n1:n2,ienergy) &
+coeff_fm2(s) * fm2(l1:l2,m1:m2,n1:n2) &
+coeff_fsub(s)*fsub(l1:l2,m1:m2,n1:n2,ienergy) &
+coeff_dfm1(s) *dt_sub_RKL* dfm1(l1:l2,m1:m2,n1:n2,ienergy) &
+coeff_dfsub(s)*dt_sub_RKL*dfsub(l1:l2,m1:m2,n1:n2,ienergy)
!
! set Yj-1 and Yj-2 for the next sub step
fm2(l1:l2,m1:m2,n1:n2)=fm1(l1:l2,m1:m2,n1:n2,ienergy)
fm1(l1:l2,m1:m2,n1:n2,ienergy)=fj(l1:l2,m1:m2,n1:n2)
endif
enddo ! end of itRKL sub steps
!
! set inital value for the next RKL step
fsub(l1:l2,m1:m2,n1:n2,ienergy)=fj(l1:l2,m1:m2,n1:n2)
!
if (notanumber(fsub(:,:,:,ienergy))) then
print*, 'fsub contains NaN in proc',iproc_world, 'in No.',j,'subcycle'
STOP
endif
!
enddo ! end of sub cycle
!
! set temperature after heat conduction back to f-array
f(l1:l2,m1:m2,n1:n2,ienergy)=fsub(l1:l2,m1:m2,n1:n2,ienergy)
!
ENDIF
!
endsubroutine time_step
!***********************************************************************
subroutine RKL_coeff(itRKL, coeff_fm1, coeff_fm2, coeff_fsub, coeff_dfm1, coeff_dfsub)
!
real, dimension (itorder) :: a,b,coeff_fm1, coeff_fm2, coeff_fsub, coeff_dfm1, coeff_dfsub
integer :: j,itRKL
!
b(1) = 1d0/3d0
b(2) = 1d0/3d0
a(1) = 1d0-b(1)
a(2) = 1d0-b(2)
!
coeff_dfm1(1) = 4d0/3d0/(itRKL*itRKL+itRKL-2d0)
!
coeff_fm1(2) = 3d0/2d0
coeff_fm2(2) = -1d0/2d0
coeff_fsub(2) = 1d0-coeff_fm1(2)-coeff_fm2(2)
coeff_dfm1(2) = 6d0/(itRKL*itRKL+itRKL-2d0)
coeff_dfsub(2) = -a(1)*coeff_dfm1(2)
!
do j=3, itRKL
!
b(j) = (j*j+j-2d0)/(2d0*j)/(j+1d0)
a(j) = 1d0-b(j)
!
coeff_fm1(j) = (2d0*j-1d0)/j*b(j)/b(j-1)
coeff_fm2(j) = -(j-1d0)/j*b(j)/b(j-2)
coeff_fsub(j) = 1d0-coeff_fm1(j)-coeff_fm2(j)
coeff_dfm1(j) = 4d0*(2d0*j-1d0)/j/(itRKL*itRKL+itRKL-2d0)*b(j)/b(j-2)
coeff_dfsub(j) = -a(j-1)*coeff_dfm1(j)
enddo
!
endsubroutine RKL_coeff
!***********************************************************************
subroutine pde_energy_only(f,df,p,dt_in)
!
! Calculate thermal conduction for sub_cycle time step.
!
! 22-jun-12/feng/bingert: coded
!
use Boundcond
use BorderProfiles, only: calc_pencils_borderprofiles
use Density
use Diagnostics
use Energy
use EquationOfState
use Grid, only: calc_pencils_grid
use Hydro
use Lorenz_gauge
use Magnetic
use Special
use Sub
!
logical :: early_finalize
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real :: dt_in
type (pencil_case) :: p
!
intent(inout) :: f ! inout due to lshift_datacube_x,
! density floor, or velocity ceiling
intent(out) :: df, p
!
! Print statements when they are first executed.
!
headtt = headt .and. lfirst .and. lroot
!
if (headtt.or.ldebug) print*,'pde: ENTER'
!
! For chemistry with LSODE
!
lchemonly=.false.
!
! Shift entire data cube by one grid point at the beginning of each
! time-step. Useful for smearing out possible x-dependent numerical
! diffusion, e.g. in a linear shear flow.
!
if (lfirst .and. lshift_datacube_x) then
call boundconds_x(f)
f=cshift(f,1,1)
endif
!
! Need to finalize communication early either for test purposes, or
! when radiation transfer of global ionization is calculated.
! This could in principle be avoided (but it not worth it now)
!
early_finalize=test_nonblocking.or. &
leos_ionization.or.lradiation_ray.or. &
lhyperviscosity_strict.or.lhyperresistivity_strict.or. &
ltestscalar.or.ltestfield.or.ltestflow.or. &
lparticles_spin.or.lsolid_cells.or. &
lchemistry.or.lweno_transport
!
! Write crash snapshots to the hard disc if the time-step is very low.
! The user must have set crash_file_dtmin_factor>0.0 in &run_pars for
! this to be done.
!
!
! For debugging purposes impose minimum or maximum value on certain variables.
!
! Sven.Bingert: routine not existing
! call impose_entropy_floor(f)
!
! Call "before_boundary" hooks (for f array precalculation)
!
if (lspecial) call special_before_boundary(f)
!
! Fetch fp to the special module
!
! Prepare x-ghost zones; required before f-array communication
! AND shock calculation
!
call boundconds_x(f,ilnTT,ilnTT)
!
! Initiate (non-blocking) communication and do boundary conditions.
! Required order:
! 1. x-boundaries (x-ghost zones will be communicated) - done above
! 2. communication
! 3. y- and z-boundaries
!
if (ldebug) print*,'pde: bef. initiate_isendrcv_bdry'
call initiate_isendrcv_bdry(f,ilnTT,ilnTT)
if (early_finalize) then
call finalize_isendrcv_bdry(f,ilnTT,ilnTT)
call boundconds_y(f,ilnTT,ilnTT)
call boundconds_z(f,ilnTT,ilnTT)
endif
!
!
!------------------------------------------------------------------------------
! Do loop over m and n.
!
mn_loop: do imn=1,ny*nz
n=nn(imn)
m=mm(imn)
lfirstpoint=(imn==1) ! true for very first m-n loop
llastpoint=(imn==(ny*nz)) ! true for very last m-n loop
!
! Make sure all ghost points are set.
!
if (.not.early_finalize.and.necessary(imn)) then
call finalize_isendrcv_bdry(f,ilnTT,ilnTT)
call boundconds_y(f,ilnTT,ilnTT)
call boundconds_z(f,ilnTT,ilnTT)
endif
!
! Grid spacing. In case of equidistant grid and cartesian coordinates
! this is calculated before the (m,n) loop.
!
if (.not. lcartesian_coords .or. .not.all(lequidist)) call get_grid_mn
!
! Calculate grid/geometry related pencils.
!
call calc_pencils_grid(f,p)
!
! Calculate pencils for the pencil_case.
! Note: some no-modules (e.g. nohydro) also calculate some pencils,
! so it would be wrong to check for lhydro etc in such cases.
!
! To check ghost cell consistency, please uncomment the following 2 lines:
! if (.not. lpencil_check_at_work .and. necessary(imn)) &
! call check_ghosts_consistency (f, 'before calc_pencils_*')
call calc_pencils_sub_cycle(f,p)
call calc_pencils_eos(f,p)
!
! --------------------------------------------------------
! NO CALLS MODIFYING PENCIL_CASE PENCILS BEYOND THIS POINT
! --------------------------------------------------------
!
! hydro, density, and entropy evolution
! Note that pressure gradient is added in denergy_dt to momentum,
! even if lentropy=.false.
!
!!! call denergy_dt(f,df,p)
if (Kpara /= 0.) call calc_heatcond_tensor(f,df,p)
if (hcond_grad /= 0.) call calc_heatcond_glnTT(df,p)
if (hcond_grad_iso /= 0.) call calc_heatcond_glnTT_iso(df,p)
if (lfilter_farray) call filter_lnTT(f,df,dt_in)
!
! End of loops over m and n.
!
headtt=.false.
enddo mn_loop
!
!
! -------------------------------------------------------------
! NO CALLS MODIFYING DF BEYOND THIS POINT (APART FROM FREEZING)
! -------------------------------------------------------------
!
! Reset lwrite_prof.
!
lwrite_prof=.false.
!
endsubroutine pde_energy_only
!***********************************************************************
subroutine calc_pencils_sub_cycle(f,p)
!
! Calculate pencil for spitzer heat conduction and
! hcond_grad and hcond_grad_iso
!
! 22 Jun F. Chen
!
use Magnetic, only: get_bext
use Messages, only: fatal_error
use Sub, only: grad, gij, curl_mn, dot2_mn, gij_etc
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
real :: B2_ext
real, dimension (nx) :: quench
logical :: luse_Bext_in_b2=.true.
real, dimension(3) :: B_ext
!
intent(inout) :: f,p
!
! Get the external magnetic field if exists.
!
if (lmagnetic) call get_bext(B_ext)
!
p%lnrho=f(l1:l2,m,n,ilnrho)
call grad(f,ilnrho,p%glnrho)
p%aa=f(l1:l2,m,n,iax:iaz)
call gij(f,iaa,p%aij,1)
call curl_mn(p%aij,p%bb,p%aa)
p%bbb=p%bb
B2_ext=B_ext(1)**2+B_ext(2)**2+B_ext(3)**2
!
if (B2_ext/=0.0) then
if (B_ext(1)/=0.0) p%bb(:,1)=p%bb(:,1)+B_ext(1)
if (B_ext(2)/=0.0) p%bb(:,2)=p%bb(:,2)+B_ext(2)
if (B_ext(3)/=0.0) p%bb(:,3)=p%bb(:,3)+B_ext(3)
endif
!
! b2 now (since 18 June 2013) includes B_ext by default.
!
if (luse_Bext_in_b2) then
if (lpencil(i_b2)) call dot2_mn(p%bb,p%b2)
else
if (lpencil(i_b2)) call dot2_mn(p%bbb,p%b2)
endif
! bunit
quench = 1.0/max(tini,sqrt(p%b2))
if (luse_Bext_in_b2) then
p%bunit(:,1) = p%bb(:,1)*quench
p%bunit(:,2) = p%bb(:,2)*quench
p%bunit(:,3) = p%bb(:,3)*quench
else
p%bunit(:,1) = p%bbb(:,1)*quench
p%bunit(:,2) = p%bbb(:,2)*quench
p%bunit(:,3) = p%bbb(:,3)*quench
endif
!
! bij, del2a, graddiva
! For non-cartesian coordinates jj is always required for del2a=graddiva-jj
!
call gij_etc(f,iaa,p%aa,p%aij,p%bij)
!
endsubroutine calc_pencils_sub_cycle
!***********************************************************************
subroutine pushpars2c(p_par)
use Syscalls, only: copy_addr
integer, parameter :: n_pars=2
integer(KIND=ikind8), dimension(n_pars) :: p_par
call copy_addr(alpha_ts,p_par(1)) ! (3)
call copy_addr(beta_ts ,p_par(2)) ! (3)
endsubroutine pushpars2c
!***********************************************************************
endmodule Timestep