!
! [Auto-generated file, so think twice before editing]
!
! A second-order timestepping module similar to RKC (Runge-Kutta-Chebyshev).
! The schemes used here are all second-order (p=2) accurate Runge-Kutta
! schemes of stage number (number of substeps) s > 2 that trade order for
! extended stability interval.
! For this file, s=10, so we have a 2nd-order, 10-step Runge-Kutta
! scheme with a critical Courant number of ~65.3 as compared to 2.513 for
! any p=s=3 Runge-Kutta scheme (like the Williamson scheme in timestep.f90).
! Here the Courant number is
! Cou = c nu dt / dx^2 ,
! where
! c = 272/45 = 6.04444
! for 6th-order central finite differences in space.
!
! This scheme uses 5N array slots (as opposed to 2N for the Williamson
! scheme in timestep.f90), irrespective of s.
! [TODO: it currently uses more, but this should be fixed...]
module Timestep
use Cparam
use Cdata
use Equ
implicit none
private
include 'timestep.h'
contains
!***********************************************************************
subroutine initialize_timestep
endsubroutine initialize_timestep
!***********************************************************************
subroutine time_step(f,df,p)
!
! Long-time-step Runge--Kutta--Chebyshev stepping, accurate to second
! order.
!
! 18-aug-08/perl: generated
!
use Mpicomm, only: mpiallreduce_max,MPI_COMM_WORLD
real, dimension (mx,my,mz,mfarray) :: f
! real, dimension (mx,my,mz,mvar) :: fn_target, fn1_target
real, dimension (mx,my,mz,mvar) :: df, dfn
type (pencil_case) :: p
real :: dt1, dt1_local
real, save :: dt1_last=0.0
double precision :: t0
integer :: iv
! Use pointers for cheaply flipping fn and fn1 after each substep
! target :: f, df
! target :: fn_target, fn1_target
! real, pointer :: f0(:,:,:,:), df0(:,:,:,:)
! real, pointer :: fn(:,:,:,:), fn1(:,:,:,:)
real, dimension(mx,my,mz,mvar) :: f0, df0, fn, fn1
! f0 => f(:,:,:,1:mvar)
f0 = f(:,:,:,1:mvar)
! fn => fn_target;
! fn1 => fn1_target;
! Step n = 0:
lfirst = .true.
t0 = t
df0 = 0.
call pde(f,df0,p)
!
! In the first time substep we need to calculate timestep dt.
! Only do this on root processor, then broadcast dt to all others.
!
if (ldt) then
dt1_local=maxval(dt1_max(1:nx))
! Timestep growth limiter
if (real(ddt) > 0.) dt1_local=max(dt1_local,dt1_last)
call mpiallreduce_max(dt1_local,dt1,MPI_COMM_WORLD)
dt=1.0/dt1
! Timestep growth limiter
if (ddt/=0.) dt1_last=dt1_local/ddt
endif
!
! IMPLEMENT ME:
! What do we need to do with dt_beta_ts?
! if (ldt) dt_beta_ts=dt*beta_ts
!
if (ip<=6) print*, 'TIMESTEP: iproc,dt=',iproc_world,dt ! same dt everywhere?
lfirst = .false.
! Step n = 1:
fn1 = f0
fn = f0 + 0.00771156039951447*dt*df0
t = t0 + 0.00771156039951447*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 2:
fn1 = -1*fn1 &
+ 2.00307692307692*fn &
+ -0.00307692307692308*f0 &
+ 0.0618824522390463*dt*dfn &
+ -0.046435603561865*dt*df0
call swap(fn, fn1)
t = t0 + 0.0308936973543626*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 3:
fn1 = -1.18094651210614*fn1 &
+ 2.365526705788*fn &
+ -0.184580193681855*f0 &
+ 0.0730798661322767*dt*dfn &
+ -0.0547538137795433*dt*df0
call swap(fn, fn1)
t = t0 + 0.0822989781283077*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 4:
fn1 = -1.23930718573423*fn1 &
+ 2.10206609605069*fn &
+ 0.137241089683541*f0 &
+ 0.0649405937902564*dt*dfn &
+ -0.0455614110277055*dt*df0
call swap(fn, fn1)
t = t0 + 0.154090293300523*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 5:
fn1 = -1.06772119115642*fn1 &
+ 2.03801206116619*fn &
+ 0.0297091299902283*f0 &
+ 0.062961727822209*dt*dfn &
+ -0.0430338340638419*dt*df0
call swap(fn, fn1)
t = t0 + 0.246093407055357*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 6:
fn1 = -1.02235335516281*fn1 &
+ 2.01274859380081*fn &
+ 0.00960476136200285*f0 &
+ 0.0621812459073009*dt*dfn &
+ -0.0418837782853428*dt*df0
call swap(fn, fn1)
t = t0 + 0.358086898262466*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 7:
fn1 = -1.00314248567622*fn1 &
+ 1.99971612021638*fn &
+ 0.00342636545984013*f0 &
+ 0.061778624612605*dt*dfn &
+ -0.0411799682884653*dt*df0
call swap(fn, fn1)
t = t0 + 0.489804047155835*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 8:
fn1 = -0.99260553114312*fn1 &
+ 1.99160679119323*fn &
+ 0.000998739949891038*f0 &
+ 0.0615280974560168*dt*dfn &
+ -0.0406510522611588*dt*df0
call swap(fn, fn1)
t = t0 + 0.64093507601914*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 9:
fn1 = -0.985705344245516*fn1 &
+ 1.9858149132024*fn &
+ -0.000109568956881862*f0 &
+ 0.0613491649302547*dt*dfn &
+ -0.0401954228616865*dt*df0
call swap(fn, fn1)
t = t0 + 0.8111297075073*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 10:
fn1 = -0.980577292130002*fn1 &
+ 1.98124561837574*fn &
+ -0.000668326245738146*f0 &
+ 0.0612080025187571*dt*dfn &
+ -0.0397688001780052*dt*df0
call swap(fn, fn1)
! Done: last fn is the updated f:
! Increase time
t = t0 + dt
! need explicit loop here, as f(:,:,:,1:mvar) = fn
! causes a `stack smashing' exception
do iv=1,mvar
f(:,:,:,iv) = fn(:,:,:,iv)
enddo
endsubroutine time_step
!***********************************************************************
subroutine swap(a, b)
!
! Swap two pointers
!
! 18-aug-08/perl: generated
!
! real, pointer :: a(:,:,:,:), b(:,:,:,:), tmp(:,:,:,:)
! tmp => a
! a => b
! b => tmp
real :: a(:,:,:,:), b(:,:,:,:), tmp(size(a,1), size(a,2), size(a,3), size(a,4))
tmp = a
a = b
b = tmp
endsubroutine swap
!***********************************************************************
subroutine pushpars2c(p_par)
use Messages, only: fatal_error
integer, parameter :: n_pars=0
integer(KIND=ikind8), dimension(:) :: p_par
call fatal_error('time_step_RKC-10','alpha_ts, beta_ts not defined')
endsubroutine pushpars2c
!***********************************************************************
endmodule Timestep