MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_dt.t
Go to the documentation of this file.
1 module mod_dt
2 
3  implicit none
4  private
5 
6 
7  public :: setdt
8 
9 contains
10 
11 
12  !>setdt - set dt for all levels between levmin and levmax.
13  !> dtpar>0 --> use fixed dtpar for all level
14  !> dtpar<=0 --> determine CFL limited timestep
15  subroutine setdt()
17  use mod_physics
18  use mod_usr_methods, only: usr_get_dt
20  use mod_comm_lib, only: mpistop
21 
22  integer :: iigrid, igrid, ncycle, ncycle2, ifile, idim
23  double precision :: dtnew, qdtnew, dtmin_mype, factor, dx^d, dxmin^d
24 
25  double precision :: dtmax, dxmin, cmax_mype
26  double precision :: a2max_mype(ndim), tco_mype, tco_global, tmax_mype, t_peak
27  double precision :: trac_alfa, trac_dmax, trac_tau, t_bott
28 
29  integer, parameter :: niter_print = 2000
30 
31  if (dtpar<=zero) then
32  dtmin_mype=bigdouble
33  cmax_mype = zero
34  a2max_mype = zero
35  tco_mype = zero
36  tmax_mype = zero
37  !$OMP PARALLEL DO PRIVATE(igrid,qdtnew,dtnew,dx^D) REDUCTION(min:dtmin_mype) REDUCTION(max:cmax_mype,a2max_mype)
38  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
39  dtnew=bigdouble
40  dx^d=rnode(rpdx^d_,igrid);
41  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
42  block=>ps(igrid)
43  if(local_timestep) then
44  ps(igrid)%dt(ixm^t)=bigdouble
45  endif
46 
47  call getdt_courant(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x,&
48  cmax_mype,a2max_mype)
49  dtnew=min(dtnew,qdtnew)
50 
51  call phys_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
52  dtnew=min(dtnew,qdtnew)
53 
54  if (associated(usr_get_dt)) then
55  call usr_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
56  dtnew = min(dtnew,qdtnew)
57  end if
58  dtmin_mype = min(dtmin_mype,dtnew)
59  end do
60  !$OMP END PARALLEL DO
61  else
62  dtmin_mype=dtpar
63  end if
64 
65  if (dtmin_mype<dtmin) then
66  write(unitterm,*)"Error: Time step too small!", dtmin_mype
67  write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
68  write(unitterm,*)"Lower limit of time step:", dtmin
69  crash=.true.
70  end if
71 
72  if (slowsteps>it-it_init+1) then
73  factor=one-(one-dble(it-it_init+1)/dble(slowsteps))**2
74  dtmin_mype=dtmin_mype*factor
75  end if
76 
77  if(final_dt_reduction)then
78  !if (dtmin_mype>time_max-global_time) then
79  ! write(unitterm,*)"WARNING final timestep artificially reduced!"
80  ! write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
81  !endif
82  if(time_max-global_time<=dtmin) then
83  !write(unitterm,*)'Forcing to leave timeloop as time is reached!'
84  final_dt_exit=.true.
85  endif
86  dtmin_mype=min(dtmin_mype,time_max-global_time)
87  end if
88 
89  if (dtpar<=zero) then
90  call mpi_allreduce(dtmin_mype,dt,1,mpi_double_precision,mpi_min, &
91  icomm,ierrmpi)
92  else
93  dt=dtmin_mype
94  end if
95 
96  if(any(dtsave(1:nfile)<bigdouble).or.any(tsave(isavet(1:nfile),1:nfile)<bigdouble))then
97  dtmax = minval(ceiling(global_time/dtsave(1:nfile))*dtsave(1:nfile))-global_time
98  do ifile=1,nfile
99  dtmax = min(tsave(isavet(ifile),ifile)-global_time,dtmax)
100  end do
101  if(dtmax > smalldouble)then
102  dt=min(dt,dtmax)
103  else
104  ! dtmax=0 means dtsave is divisible by global_time
105  dt=min(dt,minval(dtsave(1:nfile)))
106  end if
107  end if
108 
109  if(mype==0) then
110  if(any(dtsave(1:nfile)<dt)) then
111  write(unitterm,*) 'Warning: timesteps: ',dt,' exceeding output intervals ', dtsave(1:nfile)
112  endif
113  endif
114 
115  if(is_sts_initialized()) then
117  qdtnew = 0.5d0 * dt
118  if (set_dt_sts_ncycles(qdtnew)) then
119  dt = 2.d0*qdtnew
120  !a quick way to print the reduction of time only every niter_print iterations
121  !Note that niter_print is a parameter variable hardcoded to the value of 200
122  if(mype==0 .and. mod(it-1, niter_print) .eq. 0) then
123  write(*,*) 'Max number of STS cycles exceeded, reducing dt to',dt
124  endif
125  endif
126  else
127  if(set_dt_sts_ncycles(dt))then
128  if(mype==0 .and. mod(it-1, niter_print) .eq. 0) then
129  write(*,*) 'Max number of STS cycles exceeded, reducing dt to',dt
130  endif
131  endif
132  endif
133  endif
134 
135  ! global Lax-Friedrich finite difference flux splitting needs fastest wave-speed
136  ! so does GLM:
137  if(need_global_cmax) call mpi_allreduce(cmax_mype,cmax_global,1,&
138  mpi_double_precision,mpi_max,icomm,ierrmpi)
139  if(need_global_a2max) call mpi_allreduce(a2max_mype,a2max_global,ndim,&
140  mpi_double_precision,mpi_max,icomm,ierrmpi)
141 
142  ! transition region adaptive thermal conduction (Johnston 2019 ApJL, 873, L22)
143  ! transition region adaptive thermal conduction (Johnston 2020 A&A, 635, 168)
144  if(phys_trac) then
145  t_bott=2.d4/unit_temperature
146  call mpi_allreduce(tmax_mype,t_peak,1,mpi_double_precision,&
147  mpi_max,icomm,ierrmpi)
148  ! TODO trac stuff should not be here at all
149  if (phys_trac_type==1) then
150  !> 1D TRAC method
151  trac_dmax=0.1d0
152  trac_tau=1.d0/unit_time
153  trac_alfa=trac_dmax**(dt/trac_tau)
154  tco_global=zero
155  {^ifoned
156  call mpi_allreduce(tco_mype,tco_global,1,mpi_double_precision,&
157  mpi_max,icomm,ierrmpi)
158  }
159 
160  endif
161  if(.not. associated(phys_trac_after_setdt)) call mpistop("phys_trac_after_setdt not set")
162  ! trac_alfa,tco_global are set only for phys_trac_type=1, should not be a problem when not initialized
163  ! side effect of modifying T_bott from mod_trac -> T_bott sent as param
164  call phys_trac_after_setdt(tco_global,trac_alfa,t_peak, t_bott)
165 
166  end if
167 
168  contains
169 
170  !> compute CFL limited dt (for variable time stepping)
171  subroutine getdt_courant(w,ixI^L,ixO^L,dtnew,dx^D,x,cmax_mype,a2max_mype)
175 
176  integer, intent(in) :: ixI^L, ixO^L
177  double precision, intent(in) :: x(ixI^S,1:ndim)
178  double precision, intent(in) :: dx^D
179  double precision, intent(inout) :: w(ixI^S,1:nw), dtnew, cmax_mype, a2max_mype(ndim)
180 
181  integer :: idims
182  integer :: hxO^L
183  double precision :: courantmax, dxinv(1:ndim), courantmaxtot, courantmaxtots
184  double precision :: cmax(ixI^S), cmaxtot(ixI^S)
185  double precision :: a2max(ndim),tco_local,Tmax_local
186 
187  dtnew=bigdouble
188  courantmax=zero
189  courantmaxtot=zero
190  courantmaxtots=zero
191 
192  ! local timestep dt has to be calculated in the
193  ! extended region because of the calculation from the
194  ! div fluxes in mod_finite_volume
195  if(local_timestep) then
196  hxomin^d=ixomin^d-1;
197  hxomax^d=ixomax^d;
198  else
199  hxomin^d=ixomin^d;
200  hxomax^d=ixomax^d;
201  endif
202 
203  if(need_global_a2max) then
204  call phys_get_a2max(w,x,ixi^l,ixo^l,a2max)
205  do idims=1,ndim
206  a2max_mype(idims) = max(a2max_mype(idims),a2max(idims))
207  end do
208  end if
209  if(phys_trac) then
210  call phys_get_tcutoff(ixi^l,ixo^l,w,x,tco_local,tmax_local)
211  {^ifoned tco_mype=max(tco_mype,tco_local) }
212  tmax_mype=max(tmax_mype,tmax_local)
213  end if
214 
215  !if(slab_uniform) then
216  ! ^D&dxinv(^D)=one/dx^D;
217  ! do idims=1,ndim
218  ! call phys_get_cmax(w,x,ixI^L,ixO^L,idims,cmax)
219  ! if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixO^S)))
220  ! if(need_global_a2max) a2max_mype(idims) = max(a2max_mype(idims),a2max(idims))
221  ! cmaxtot(ixO^S)=cmaxtot(ixO^S)+cmax(ixO^S)*dxinv(idims)
222  ! courantmax=max(courantmax,maxval(cmax(ixO^S)*dxinv(idims)))
223  ! courantmaxtot=courantmaxtot+courantmax
224  ! end do
225  !else
226  ! do idims=1,ndim
227  ! call phys_get_cmax(w,x,ixI^L,ixO^L,idims,cmax)
228  ! if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixO^S)))
229  ! if(need_global_a2max) a2max_mype(idims) = max(a2max_mype(idims),a2max(idims))
230  ! tmp(ixO^S)=cmax(ixO^S)/block%ds(ixO^S,idims)
231  ! cmaxtot(ixO^S)=cmaxtot(ixO^S)+tmp(ixO^S)
232  ! courantmax=max(courantmax,maxval(tmp(ixO^S)))
233  ! courantmaxtot=courantmaxtot+courantmax
234  ! end do
235  !end if
236 
237  ! these are also calculated in hxO because of local timestep
238  if(nwaux>0) call phys_get_auxiliary(ixi^l,hxo^l,w,x)
239 
240  select case (type_courant)
241  case (type_maxsum)
242  cmaxtot(hxo^s)=zero
243  if(slab_uniform) then
244  ^d&dxinv(^d)=one/dx^d;
245  do idims=1,ndim
246  call phys_get_cmax(w,x,ixi^l,hxo^l,idims,cmax)
247  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
248  cmaxtot(hxo^s)=cmaxtot(hxo^s)+cmax(hxo^s)*dxinv(idims)
249  end do
250  else
251  do idims=1,ndim
252  call phys_get_cmax(w,x,ixi^l,hxo^l,idims,cmax)
253  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
254  cmaxtot(hxo^s)=cmaxtot(hxo^s)+cmax(hxo^s)/block%ds(hxo^s,idims)
255  end do
256  end if
257  ! courantmaxtots='max(summed c/dx)'
258  courantmaxtots=max(courantmaxtots,maxval(cmaxtot(ixo^s)))
259  if (courantmaxtots>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtots)
260  if(local_timestep) then
261  block%dt(hxo^s) = courantpar/cmaxtot(hxo^s)
262  endif
263 
264  case (type_summax)
265  !TODO this should be mod_input_output?
266  if(local_timestep) then
267  call mpistop("Type courant summax incompatible with local_timestep")
268  endif
269  if(slab_uniform) then
270  ^d&dxinv(^d)=one/dx^d;
271  do idims=1,ndim
272  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
273  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
274  courantmax=max(courantmax,maxval(cmax(ixo^s)*dxinv(idims)))
275  courantmaxtot=courantmaxtot+courantmax
276  end do
277  else
278  do idims=1,ndim
279  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
280  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
281  courantmax=max(courantmax,maxval(cmax(ixo^s)/block%ds(ixo^s,idims)))
282  courantmaxtot=courantmaxtot+courantmax
283  end do
284  end if
285  ! courantmaxtot='summed max(c/dx)'
286  if (courantmaxtot>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtot)
287  case (type_minimum)
288  if(local_timestep) then
289  call mpistop("Type courant not implemented for local_timestep, use maxsum")
290  endif
291  if(slab_uniform) then
292  ^d&dxinv(^d)=one/dx^d;
293  do idims=1,ndim
294  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
295  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
296  courantmax=max(courantmax,maxval(cmax(ixo^s)*dxinv(idims)))
297  end do
298  else
299  do idims=1,ndim
300  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
301  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
302  courantmax=max(courantmax,maxval(cmax(ixo^s)/block%ds(ixo^s,idims)))
303  end do
304  end if
305  ! courantmax='max(c/dx)'
306  if (courantmax>smalldouble) dtnew=min(dtnew,courantpar/courantmax)
307  end select
308 
309  end subroutine getdt_courant
310 
311  end subroutine setdt
312 
313 end module mod_dt
subroutine getdt_courant(w, ixIL, ixOL, dtnew, dxD, x, cmax_mype, a2max_mype)
compute CFL limited dt (for variable time stepping)
Definition: mod_dt.t:172
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Definition: mod_dt.t:1
subroutine, public setdt()
setdt - set dt for all levels between levmin and levmax. dtpar>0 --> use fixed dtpar for all level dt...
Definition: mod_dt.t:16
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision unit_time
Physical scaling factor for time.
double precision global_time
The global simulation time.
double precision time_max
End time for the simulation.
integer it
Number of time steps taken.
integer it_init
initial iteration count
integer, parameter type_maxsum
integer switchers for type courant
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
double precision dtpar
If dtpar is positive, it sets the timestep dt, otherwise courantpar is used to limit the time step ba...
integer, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
logical need_global_a2max
global value for schmid scheme
double precision courantpar
The Courant (CFL) number used for the simulation.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
integer slowsteps
If > 1, then in the first slowsteps-1 time steps dt is reduced by a factor .
integer type_courant
How to compute the CFL-limited time step.
integer, parameter unitterm
Unit for standard output.
double precision, dimension(nfile) dtsave
Repeatedly save output of type N when dtsave(N) simulation time has passed.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
logical final_dt_reduction
If true, allow final dt reduction for matching time_max on output.
integer, parameter type_summax
double precision, dimension(:,:), allocatable dx
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
double precision, dimension(nsavehi, nfile) tsave
Save output of type N on times tsave(:, N)
logical need_global_cmax
need global maximal wave speed
logical crash
Save a snapshot before crash a run met unphysical values.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
integer, parameter nfile
Number of output methods.
logical final_dt_exit
Force timeloop exit when final dt < dtmin.
integer, parameter type_minimum
double precision, dimension(ndim) dxlevel
integer, dimension(nfile) isavet
double precision, dimension(ndim) a2max_global
global largest a2 for schmid scheme
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_get_a2max), pointer phys_get_a2max
Definition: mod_physics.t:59
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:66
procedure(sub_get_tcutoff), pointer phys_get_tcutoff
Definition: mod_physics.t:60
procedure(sub_get_auxiliary), pointer phys_get_auxiliary
Definition: mod_physics.t:90
procedure(sub_trac_after_setdt), pointer phys_trac_after_setdt
Definition: mod_physics.t:61
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:58
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
integer, public sourcetype_sts
pure logical function, public is_sts_initialized()
logical function, public set_dt_sts_ncycles(my_dt)
This sets the explicit dt and calculates the number of cycles for each of the terms implemented with ...
integer, parameter, public sourcetype_sts_split
Module with all the methods that users can customize in AMRVAC.
procedure(get_dt), pointer usr_get_dt