!
! [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=25, so we have a 2nd-order, 25-step Runge-Kutta
! scheme with a critical Courant number of ~408.125 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 Cdata
use Equ
implicit none
private
include 'timestep.h'
contains
!***********************************************************************
subroutine initialize_timestep
if (dt0 < 0.) dt = 0
ldt = (dt==0.)
lcourant_dt = .true.
num_substeps = 25
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
real(KIND=rkind8) :: 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.00122590002096801*dt*df0
t = t0 + 0.00122590002096801*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 2:
fn1 = -1*fn1 &
+ 2.00049230769231*fn &
+ -0.000492307692307692*f0 &
+ 0.00981202892206177*dt*dfn &
+ -0.00735962536011542*dt*df0
call swap(fn, fn1)
t = t0 + 0.0049048071238927*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 3:
fn1 = -1.18450487013592*fn1 &
+ 2.36959288113099*fn &
+ -0.185088010995069*f0 &
+ 0.0116223960440967*dt*dfn &
+ -0.00871536711232105*dt*df0
call swap(fn, fn1)
t = t0 + 0.013077340430736*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 4:
fn1 = -1.24827889022172*fn1 &
+ 2.10819928284186*fn &
+ 0.140079607379862*f0 &
+ 0.0103403108610682*dt*dfn &
+ -0.00727299814858042*dt*df0
call swap(fn, fn1)
t = t0 + 0.0245143852366801*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 5:
fn1 = -1.07801830248992*fn1 &
+ 2.04638924467856*fn &
+ 0.0316290578113557*f0 &
+ 0.0100371445455568*dt*dfn &
+ -0.00689405919497779*dt*df0
call swap(fn, fn1)
t = t0 + 0.0392114473066032*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 6:
fn1 = -1.03465987422963*fn1 &
+ 2.02340646703125*fn &
+ 0.0112534071983794*f0 &
+ 0.00992441845402541*dt*dfn &
+ -0.00673830461455823*dt*df0
call swap(fn, fn1)
t = t0 + 0.0571627591555031*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 7:
fn1 = -1.01760320970258*fn1 &
+ 2.01265393155044*fn &
+ 0.00494927815213762*f0 &
+ 0.00987167934139925*dt*dfn &
+ -0.00665754853714671*dt*df0
call swap(fn, fn1)
t = t0 + 0.0783612881205638*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 8:
fn1 = -1.00925537641767*fn1 &
+ 2.00680760050515*fn &
+ 0.00244777591252246*f0 &
+ 0.00984300421523965*dt*dfn &
+ -0.00660846461060895*dt*df0
call swap(fn, fn1)
t = t0 + 0.102798746181792*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 9:
fn1 = -1.0045410436821*fn1 &
+ 2.0032526238769*fn &
+ 0.00128841980519865*f0 &
+ 0.00982556774054818*dt*dfn &
+ -0.0065746941031726*dt*df0
call swap(fn, fn1)
t = t0 + 0.130465601504401*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 10:
fn1 = -1.00157786203203*fn1 &
+ 2.00088794645573*fn &
+ 0.000689915576302148*f0 &
+ 0.0098139694538871*dt*dfn &
+ -0.00654897883927557*dt*df0
call swap(fn, fn1)
t = t0 + 0.161351091671532*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 11:
fn1 = -0.999548016840505*fn1 &
+ 1.99919273711209*fn &
+ 0.000355279728414495*f0 &
+ 0.00980565477902192*dt*dfn &
+ -0.00652771355278923*dt*df0
call swap(fn, fn1)
t = t0 + 0.195443238571438*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 12:
fn1 = -0.998054466378356*fn1 &
+ 1.99789816800909*fn &
+ 0.000156298369269838*f0 &
+ 0.00979930516726314*dt*dfn &
+ -0.00650895086212785*dt*df0
call swap(fn, fn1)
t = t0 + 0.232728864898978*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 13:
fn1 = -0.996887746326627*fn1 &
+ 1.9968556920587*fn &
+ 3.20542679236489e-05*f0 &
+ 0.00979419202379521*dt*dfn &
+ -0.00649156874230619*dt*df0
call swap(fn, fn1)
t = t0 + 0.273193612227117*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 14:
fn1 = -0.995929939812475*fn1 &
+ 1.99597859392767*fn &
+ -4.86541151938653e-05*f0 &
+ 0.00978989002663376*dt*dfn &
+ -0.00647488701804963*dt*df0
call swap(fn, fn1)
t = t0 + 0.31682196060022*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 15:
fn1 = -0.995111103778451*fn1 &
+ 1.99521391252177*fn &
+ -0.000102808743319223*f0 &
+ 0.0097861394118266*dt*dfn &
+ -0.00645847638805497*dt*df0
call swap(fn, fn1)
t = t0 + 0.363597249597207*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 16:
fn1 = -0.99438792546303*fn1 &
+ 1.99452805395108*fn &
+ -0.000140128488054702*f0 &
+ 0.00978277540782312*dt*dfn &
+ -0.00644205711742458*dt*df0
call swap(fn, fn1)
t = t0 + 0.4135017008091*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 17:
fn1 = -0.993732597180402*fn1 &
+ 1.99389901310145*fn &
+ -0.000166415921044833*f0 &
+ 0.0097796900837826*dt*dfn &
+ -0.00642544240066078*dt*df0
call swap(fn, fn1)
t = t0 + 0.466516441672272*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 18:
fn1 = -0.993126709747659*fn1 &
+ 1.9933119727985*fn &
+ -0.0001852630508383*f0 &
+ 0.00977681076432268*dt*dfn &
+ -0.00640850526176165*dt*df0
call swap(fn, fn1)
t = t0 + 0.522621530595651*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 19:
fn1 = -0.992557750743432*fn1 &
+ 1.99275671495773*fn &
+ -0.000198964214295297*f0 &
+ 0.00977408733170967*dt*dfn &
+ -0.00639115845924789*dt*df0
call swap(fn, fn1)
t = t0 + 0.581795983317371*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 20:
fn1 = -0.992017019612368*fn1 &
+ 1.99222604521715*fn &
+ -0.000209025604779197*f0 &
+ 0.00977148449898564*dt*dfn &
+ -0.00637334187594473*dt*df0
call swap(fn, fn1)
t = t0 + 0.644017800423876*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 21:
fn1 = -0.991498345148503*fn1 &
+ 1.99171480539336*fn &
+ -0.000216460244862063*f0 &
+ 0.00976897696625593*dt*dfn &
+ -0.00635501437190806*dt*df0
call swap(fn, fn1)
t = t0 + 0.709263995962232*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 22:
fn1 = -0.990997273618344*fn1 &
+ 1.99121923793624*fn &
+ -0.000221964317900269*f0 &
+ 0.00976654630346188*dt*dfn &
+ -0.00633614838280148*dt*df0
call swap(fn, fn1)
t = t0 + 0.777510627074447*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 23:
fn1 = -0.990510541600783*fn1 &
+ 1.99073656721709*fn &
+ -0.000226025616309056*f0 &
+ 0.00976417889668009*dt*dfn &
+ -0.00631672625424829*dt*df0
call swap(fn, fn1)
t = t0 + 0.84873282458094*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 24:
fn1 = -0.990035725853776*fn1 &
+ 1.99026471780727*fn &
+ -0.00022899195349787*f0 &
+ 0.00976186456633341*dt*dfn &
+ -0.00629673770113825*dt*df0
call swap(fn, fn1)
t = t0 + 0.922904824438875*dt
f(:,:,:,1:mvar) = fn
dfn = 0.
call pde(f,dfn,p)
! Step n = 25:
fn1 = -0.989571005978995*fn1 &
+ 1.9898021212935*fn &
+ -0.000231115314500454*f0 &
+ 0.00975959561965714*dt*dfn &
+ -0.0062761780121769*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-25','alpha_ts, beta_ts not defined')
endsubroutine pushpars2c
!***********************************************************************
endmodule Timestep