!
! [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=20, so we have a 2nd-order, 20-step Runge-Kutta
! scheme with a critical Courant number of ~261.2 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.00191678900002874*dt*df0
t = t0 + 0.00191678900002874*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 2:
fn1 = -1*fn1 &
+ 2.00076923076923*fn &
+ -0.000769230769230769*f0 &
+ 0.0153461098932349*dt*dfn &
+ -0.0115110574401004*dt*df0
call swap(fn, fn1)
t = t0 + 0.00767010490626887*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 3:
fn1 = -1.18412255004256*fn1 &
+ 2.36915596358514*fn &
+ -0.185033413542588*f0 &
+ 0.0181716747800106*dt*dfn &
+ -0.0136252635472204*dt*df0
call swap(fn, fn1)
t = t0 + 0.0204483729341932*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 4:
fn1 = -1.2473124920472*fn1 &
+ 2.10753900020943*fn &
+ 0.13977349183777*f0 &
+ 0.0161650452256593*dt*dfn &
+ -0.0113668150740147*dt*df0
call swap(fn, fn1)
t = t0 + 0.038326955936382*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 5:
fn1 = -1.07690647956537*fn1 &
+ 2.04548533310867*fn &
+ 0.0314211464566949*f0 &
+ 0.0156890870891778*dt*dfn &
+ -0.0107704374044694*dt*df0
call swap(fn, fn1)
t = t0 + 0.0612948906058786*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 6:
fn1 = -1.03332741777686*fn1 &
+ 2.02225344008458*fn &
+ 0.0110739776922838*f0 &
+ 0.0155108960325118*dt*dfn &
+ -0.0105222825595372*dt*df0
call swap(fn, fn1)
t = t0 + 0.0893381224513312*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 7:
fn1 = -1.01603262651939*fn1 &
+ 2.01125006435675*fn &
+ 0.00478256216264779*f0 &
+ 0.0154264989863563*dt*dfn &
+ -0.0103906583758684*dt*df0
call swap(fn, fn1)
t = t0 + 0.122439536445729*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 8:
fn1 = -1.00744075365396*fn1 &
+ 2.00515266595057*fn &
+ 0.00228808770339126*f0 &
+ 0.0153797313009256*dt*dfn &
+ -0.0103078345924054*dt*df0
call swap(fn, fn1)
t = t0 + 0.160578994218027*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 9:
fn1 = -1.00248035919761*fn1 &
+ 2.00134715579027*fn &
+ 0.0011332034073402*f0 &
+ 0.0153505426387742*dt*dfn &
+ -0.0102482478907792*dt*df0
call swap(fn, fn1)
t = t0 + 0.203733377629801*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 10:
fn1 = -0.999270803493343*fn1 &
+ 1.99873293999996*fn &
+ 0.000537863493382398*f0 &
+ 0.0153304913294149*dt*dfn &
+ -0.0102005650461667*dt*df0
call swap(fn, fn1)
t = t0 + 0.251876638552983*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 11:
fn1 = -0.99699517423012*fn1 &
+ 1.99678950821599*fn &
+ 0.000205666014129042*f0 &
+ 0.0153155850037537*dt*dfn &
+ -0.0101591652646928*dt*df0
call swap(fn, fn1)
t = t0 + 0.304979854639854*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 12:
fn1 = -0.995256978396903*fn1 &
+ 1.99524827993558*fn &
+ 8.69846132574327e-06*f0 &
+ 0.0153037636211584*dt*dfn &
+ -0.0101210209382359*dt*df0
call swap(fn, fn1)
t = t0 + 0.363011290853083*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 13:
fn1 = -0.993847125881325*fn1 &
+ 1.99396091341852*fn &
+ -0.000113787537197651*f0 &
+ 0.0152938893849198*dt*dfn &
+ -0.0100843959883698*dt*df0
call swap(fn, fn1)
t = t0 + 0.425936466501677*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 14:
fn1 = -0.992647979448224*fn1 &
+ 1.99284087260328*fn &
+ -0.000192893155050772*f0 &
+ 0.0152852985543675*dt*dfn &
+ -0.010048246240369*dt*df0
call swap(fn, fn1)
t = t0 + 0.493718227508551*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 15:
fn1 = -0.991589820209969*fn1 &
+ 1.99183535813111*fn &
+ -0.000245537921142193*f0 &
+ 0.0152775861528812*dt*dfn &
+ -0.0100119208019466*dt*df0
call swap(fn, fn1)
t = t0 + 0.56631682361702*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 16:
fn1 = -0.990629524626519*fn1 &
+ 1.99091092510939*fn &
+ -0.000281400482868114*f0 &
+ 0.0152704956546258*dt*dfn &
+ -0.00997500362374063*dt*df0
call swap(fn, fn1)
t = t0 + 0.643689990227013*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 17:
fn1 = -0.989739450669471*fn1 &
+ 1.99004570747471*fn &
+ -0.00030625680523677*f0 &
+ 0.0152638593446011*dt*dfn &
+ -0.00993722489141483*dt*df0
call swap(fn, fn1)
t = t0 + 0.725793034537301*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 18:
fn1 = -0.988901337853992*fn1 &
+ 1.98922501882397*fn &
+ -0.000323680969975861*f0 &
+ 0.0152575645765545*dt*dfn &
+ -0.00989840921008424*dt*df0
call swap(fn, fn1)
t = t0 + 0.812578925657521*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 19:
fn1 = -0.988102809471123*fn1 &
+ 1.98843876457778*fn &
+ -0.000335955106657913*f0 &
+ 0.0152515339240032*dt*dfn &
+ -0.00985844411613771*dt*df0
call swap(fn, fn1)
t = t0 + 0.903998388343368*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 20:
fn1 = -0.987335290155306*fn1 &
+ 1.98767986723567*fn &
+ -0.000344577080366816*f0 &
+ 0.0152457131017762*dt*dfn &
+ -0.00981726028722926*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-20','alpha_ts, beta_ts not defined')
endsubroutine pushpars2c
!***********************************************************************
endmodule Timestep