MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
setdt.t
Go to the documentation of this file.
1 !>setdt - set dt for all levels between levmin and levmax.
2 !> dtpar>0 --> use fixed dtpar for all level
3 !> dtpar<=0 --> determine CFL limited timestep
4 subroutine setdt()
6  use mod_physics
7  use mod_usr_methods, only: usr_get_dt
9 
10  integer :: iigrid, igrid, ncycle, ncycle2, ifile, idim
11  double precision :: dtnew, qdtnew, dtmin_mype, factor, dx^D, dxmin^D
12 
13  double precision :: dtmax, dxmin, cmax_mype, v(ixg^t)
14 
15  if (dtpar<=zero) then
16  dtmin_mype=bigdouble
17  cmax_mype = zero
18  !$OMP PARALLEL DO PRIVATE(igrid,qdtnew,dtnew,dx^D)
19  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
20  dtnew=bigdouble
21  dx^d=rnode(rpdx^d_,igrid);
22  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
23  saveigrid = igrid
24  block=>ps(igrid)
25  block%iw0=0
26 
27  if (nwaux>0) then
28  call phys_get_aux(.true.,ps(igrid)%w,&
29  ps(igrid)%x,ixg^ll,ixm^ll,'setdt')
30  end if
31 
32  call getdt_courant(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,ps(igrid)%x)
33  dtnew=min(dtnew,qdtnew)
34 
35  call phys_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
36  dtnew=min(dtnew,qdtnew)
37 
38  if (associated(usr_get_dt)) then
39  call usr_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
40  end if
41 
42  dtnew = min(dtnew,qdtnew)
43  dtmin_mype = min(dtmin_mype,dtnew)
44  dt_grid(igrid) = dtnew
45  end do
46  !$OMP END PARALLEL DO
47  else
48  dtmin_mype=dtpar
49  end if
50 
51  if (dtmin_mype<dtmin) then
52  write(unitterm,*)"Error: Time step too small!", dtmin_mype
53  write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
54  write(unitterm,*)"Lower limit of time step:", dtmin
55  crash=.true.
56  end if
57 
58  if (slowsteps>it-it_init+1) then
59  factor=one-(one-dble(it-it_init+1)/dble(slowsteps))**2
60  dtmin_mype=dtmin_mype*factor
61  end if
62 
63 
64  dtmin_mype=min(dtmin_mype,time_max-global_time)
65 
66  if (dtpar<=zero) then
67  call mpi_allreduce(dtmin_mype,dt,1,mpi_double_precision,mpi_min, &
68  icomm,ierrmpi)
69  else
70  dt=dtmin_mype
71  end if
72 
73  if(any(dtsave(1:nfile)<bigdouble).or.any(tsave(isavet(1:nfile),1:nfile)<bigdouble))then
74  dtmax = minval(ceiling(global_time/dtsave(1:nfile))*dtsave(1:nfile))-global_time
75  do ifile=1,nfile
76  dtmax = min(tsave(isavet(ifile),ifile)-global_time,dtmax)
77  end do
78  if(dtmax > smalldouble)then
79  dt=min(dt,dtmax)
80  else
81  ! dtmax=0 means dtsave is divisible by global_time
82  dt=min(dt,minval(dtsave(1:nfile)))
83  end if
84  end if
85 
86  if(mype==0) then
87  if(any(dtsave(1:nfile)<dt)) then
88  write(unitterm,*) 'Warning: timesteps: ',dt,' exceeding output intervals ', dtsave(1:nfile)
89  endif
90  endif
91 
92  ! estimate time step of thermal conduction
93  if(associated(phys_getdt_heatconduct)) then
94  dtmin_mype=bigdouble
95  !$OMP PARALLEL DO PRIVATE(igrid,qdtnew,&
96  !$OMP& dx^D) REDUCTION(min:dtmin_mype)
97  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
98  dx^d=rnode(rpdx^d_,igrid);
99  saveigrid = igrid
100  block=>ps(igrid)
101  qdtnew=bigdouble
102  call phys_getdt_heatconduct(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
103  dtmin_mype=min(dtmin_mype,qdtnew)
104  end do
105  !$OMP END PARALLEL DO
106  call mpi_allreduce(dtmin_mype,dtnew,1,mpi_double_precision,mpi_min, &
107  icomm,ierrmpi)
108  if(all(flux_scheme=='nul')) dt=min(dt,dtnew)
109  ncycle=ceiling(dt/dtnew)
110  if (ncycle>tc_ncycles) then
111  if(mype==0 .and. .false.) then
112  write(*,*) 'CLF time step is too many times larger than conduction time step',ncycle
113  write(*,*) 'reducing dt to',tc_ncycles,'times of dt_impl!!'
114  endif
115  dt=tc_ncycles*dtnew
116  endif
117  ! get number of sub-steps of supertime stepping (Meyer 2012 MNRAS 422,2102)
118  if(dt/dtnew< 0.5d0) then
119  s=1
120  else if(dt/dtnew< 2.d0) then
121  s=2
122  else
123  s=ceiling((dsqrt(9.d0+8.d0*dt/dtnew)-1.d0)/2.d0)
124  ! only use odd s number
125  s=s/2*2+1
126  endif
127  dt_tc=dt*0.5d0
128  if(mype==0 .and. .false.) write(*,*) 'supertime steps:',s,' normal subcycles:',&
129  ceiling(dt/dtnew/2.d0)
130  endif
131 
132  !$OMP PARALLEL DO PRIVATE(igrid)
133  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
134  dt_grid(igrid)=dt
135  end do
136  !$OMP END PARALLEL DO
137 
138 
139  ! global Lax-Friedrich finite difference flux splitting needs fastest wave-speed
140  ! so does GLM:
141  if(need_global_cmax) call mpi_allreduce(cmax_mype,cmax_global,1,&
142  mpi_double_precision,mpi_max,icomm,ierrmpi)
143 
144  contains
145 
146  !> compute CFL limited dt (for variable time stepping)
147  subroutine getdt_courant(w,ixI^L,ixO^L,dtnew,x)
149  use mod_physics, only: phys_get_cmax
150 
151  integer, intent(in) :: ixI^L, ixO^L
152  double precision, intent(in) :: x(ixi^s,1:ndim)
153  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
154 
155  integer :: idims
156  double precision :: courantmax, dxinv(1:ndim), courantmaxtot, courantmaxtots
157  double precision :: cmax(ixi^s), cmaxtot(ixi^s), tmp(ixi^s)
158 
159  dtnew=bigdouble
160 
161  courantmax=zero
162  courantmaxtot=zero
163  courantmaxtots=zero
164 
165  ^d&dxinv(^d)=one/dx^d;
166 
167  cmaxtot(ixo^s)=zero
168 
169  do idims=1,ndim
170  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
171  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
172  if(slab_uniform) then
173  cmaxtot(ixo^s)=cmaxtot(ixo^s)+cmax(ixo^s)*dxinv(idims)
174  courantmax=max(courantmax,maxval(cmax(ixo^s)*dxinv(idims)))
175  else
176  tmp(ixo^s)=cmax(ixo^s)/block%ds(ixo^s,idims)
177  cmaxtot(ixo^s)=cmaxtot(ixo^s)+tmp(ixo^s)
178  courantmax=max(courantmax,maxval(tmp(ixo^s)))
179  end if
180  courantmaxtot=courantmaxtot+courantmax
181  end do
182 
183  select case (typecourant)
184  case ('minimum')
185  ! courantmax='max(c/dx)'
186  if (courantmax>smalldouble) dtnew=min(dtnew,courantpar/courantmax)
187  case ('summax')
188  ! courantmaxtot='summed max(c/dx)'
189  if (courantmaxtot>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtot)
190  case ('maxsum')
191  ! courantmaxtots='max(summed c/dx)'
192  courantmaxtots=max(courantmaxtots,maxval(cmaxtot(ixo^s)))
193  if (courantmaxtots>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtots)
194  case default
195  write(unitterm,*)'Unknown typecourant=',typecourant
196  call mpistop("Error from getdt_courant: no such typecourant!")
197  end select
198 
199  end subroutine getdt_courant
200 
201 end subroutine setdt
This module contains definitions of global parameters and variables and some generic functions/subrou...
procedure(get_dt), pointer usr_get_dt
double precision time_max
End time for the simulation.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
logical need_global_cmax
need global maximal wave speed
integer, public s
Number of sub-steps of supertime stepping.
integer it_init
initial iteration count
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
double precision courantpar
The Courant (CFL) number used for the simulation.
character(len=std_len), dimension(:), allocatable flux_scheme
Which spatial discretization to use (per grid level)
double precision dt_tc
Time step of thermal conduction.
integer, dimension(nfile) isavet
Module with all the methods that users can customize in AMRVAC.
procedure(getdt_heatconduct), pointer phys_getdt_heatconduct
integer tc_ncycles
Maximal substeps of TC within one fluid time step to limit fluid time step.
Thermal conduction for HD and MHD.
integer ierrmpi
A global MPI error return code.
integer ixm
the mesh range (within a block with ghost cells)
double precision dtpar
If dtpar is positive, it sets the timestep dt, otherwise courantpar is used to limit the time step ba...
double precision, dimension(:), allocatable dt_grid
integer, parameter unitterm
Unit for standard output.
integer mype
The rank of the current MPI task.
double precision global_time
The global simulation time.
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:53
double precision, dimension(nfile) dtsave
Repeatedly save output of type N when dtsave(N) simulation time has passed.
subroutine getdt_courant(w, ixIL, ixOL, dtnew, x)
compute CFL limited dt (for variable time stepping)
Definition: setdt.t:148
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
double precision, dimension(ndim) dxlevel
double precision, dimension(nsavehi, nfile) tsave
Save output of type N on times tsave(:, N)
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical crash
Save a snapshot before crash a run met unphysical values.
integer slowsteps
If > 1, then in the first slowsteps-1 time steps dt is reduced by a factor .
integer icomm
The MPI communicator.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:49
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
integer, parameter nfile
Number of output methods.
subroutine setdt()
setdt - set dt for all levels between levmin and levmax. dtpar>0 –> use fixed dtpar for all level dt...
Definition: setdt.t:5
procedure(sub_get_aux), pointer phys_get_aux
Definition: mod_physics.t:57
character(len=std_len) typecourant
How to compute the CFL-limited time step.
integer it
Number of time steps taken.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.