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  implicit none
3  private
4  public :: setdt
5 
6 contains
7  !>setdt - set dt for all levels between levmin and levmax.
8  !> dtpar>0 --> use fixed dtpar for all level
9  !> dtpar<=0 --> determine CFL limited timestep
10  subroutine setdt()
12  use mod_physics
13  use mod_usr_methods, only: usr_get_dt
15  use mod_comm_lib, only: mpistop
16 
17  integer :: iigrid, igrid, ncycle, ncycle2, ifile, idim
18  double precision :: dtnew, qdtnew, dtmin_mype, factor, dx^d, dxmin^d
19  double precision :: dtmax, dxmin, cmax_mype
20  double precision :: a2max_mype(ndim), cs2max_mype, tco_mype, tco_global, tmax_mype, t_peak
21  double precision :: trac_alfa, trac_dmax, trac_tau, t_bott
22  integer, parameter :: niter_print = 2000
23 
24  if (dtpar<=zero) then
25  dtmin_mype = bigdouble
26  cmax_mype = zero
27  a2max_mype = zero
28  cs2max_mype = zero
29  tco_mype = zero
30  tmax_mype = zero
31  !$OMP PARALLEL DO PRIVATE(igrid,qdtnew,dtnew,dx^D) REDUCTION(min:dtmin_mype) REDUCTION(max:cmax_mype,a2max_mype,cs2max_mype)
32  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
33  dtnew=bigdouble
34  dx^d=rnode(rpdx^d_,igrid);
35  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
36  block=>ps(igrid)
37  if(local_timestep) then
38  ps(igrid)%dt(ixm^t)=bigdouble
39  endif
40  call getdt_courant(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x,&
41  cmax_mype,a2max_mype,cs2max_mype)
42  dtnew=min(dtnew,qdtnew)
43 
44  call phys_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
45  dtnew=min(dtnew,qdtnew)
46 
47  if (associated(usr_get_dt)) then
48  call usr_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
49  dtnew = min(dtnew,qdtnew)
50  end if
51  dtmin_mype = min(dtmin_mype,dtnew)
52  end do
53  !$OMP END PARALLEL DO
54  else
55  dtmin_mype=dtpar
56  end if
57 
58  if (dtmin_mype<dtmin) then
59  write(unitterm,*)"Error: Time step too small!", dtmin_mype
60  write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
61  write(unitterm,*)"Lower limit of time step:", dtmin
62  crash=.true.
63  end if
64 
65  if (slowsteps>it-it_init+1) then
66  factor=one-(one-dble(it-it_init+1)/dble(slowsteps))**2
67  dtmin_mype=dtmin_mype*factor
68  end if
69 
70  if(final_dt_reduction)then
71  !if (dtmin_mype>time_max-global_time) then
72  ! write(unitterm,*)"WARNING final timestep artificially reduced!"
73  ! write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
74  !endif
75  if(time_max-global_time<=dtmin) then
76  !write(unitterm,*)'Forcing to leave timeloop as time is reached!'
77  final_dt_exit=.true.
78  endif
79  dtmin_mype=min(dtmin_mype,time_max-global_time)
80  end if
81 
82  if (dtpar<=zero) then
83  call mpi_allreduce(dtmin_mype,dt,1,mpi_double_precision,mpi_min, &
84  icomm,ierrmpi)
85  else
86  dt=dtmin_mype
87  end if
88 
89  if(any(dtsave(1:nfile)<bigdouble).or.any(tsave(isavet(1:nfile),1:nfile)<bigdouble))then
90  dtmax = minval(ceiling(global_time/dtsave(1:nfile))*dtsave(1:nfile))-global_time
91  do ifile=1,nfile
92  dtmax = min(tsave(isavet(ifile),ifile)-global_time,dtmax)
93  end do
94  if(dtmax > smalldouble)then
95  dt=min(dt,dtmax)
96  else
97  ! dtmax=0 means dtsave is divisible by global_time
98  dt=min(dt,minval(dtsave(1:nfile)))
99  end if
100  end if
101 
102  if(mype==0) then
103  if(any(dtsave(1:nfile)<dt)) then
104  write(unitterm,*) 'Warning: timesteps: ',dt,' exceeding output intervals ', dtsave(1:nfile)
105  endif
106  endif
107  if(is_sts_initialized()) then
109  qdtnew = 0.5d0 * dt
110  if (set_dt_sts_ncycles(qdtnew)) then
111  dt = 2.d0*qdtnew
112  !a quick way to print the reduction of time only every niter_print iterations
113  !Note that niter_print is a parameter variable hardcoded to the value of 200
114  if(mype==0 .and. mod(it-1, niter_print) .eq. 0) then
115  write(*,*) 'Max number of STS cycles exceeded, reducing dt to',dt
116  endif
117  endif
118  else
119  if(set_dt_sts_ncycles(dt))then
120  if(mype==0 .and. mod(it-1, niter_print) .eq. 0) then
121  write(*,*) 'Max number of STS cycles exceeded, reducing dt to',dt
122  endif
123  endif
124  endif
125  endif
126 
127  ! global Lax-Friedrich finite difference flux splitting needs fastest wave-speed
128  ! so does GLM:
129  if(need_global_cmax) call mpi_allreduce(cmax_mype, cmax_global, 1,&
130  mpi_double_precision,mpi_max,icomm,ierrmpi)
131  if(need_global_a2max) call mpi_allreduce(a2max_mype, a2max_global, ndim,&
132  mpi_double_precision,mpi_max,icomm,ierrmpi)
133  if(need_global_cs2max) call mpi_allreduce(cs2max_mype, cs2max_global, 1,&
134  mpi_double_precision,mpi_max,icomm,ierrmpi)
135 
136  ! transition region adaptive thermal conduction (Johnston 2019 ApJL, 873, L22)
137  ! transition region adaptive thermal conduction (Johnston 2020 A&A, 635, 168)
138  if(phys_trac) then
139  t_bott=2.d4/unit_temperature
140  call mpi_allreduce(tmax_mype,t_peak,1,mpi_double_precision,&
141  mpi_max,icomm,ierrmpi)
142  ! TODO trac stuff should not be here at all
143  if(phys_trac_type==1) then
144  !> 1D TRAC method
145  trac_dmax=0.1d0
146  trac_tau=1.d0/unit_time
147  trac_alfa=trac_dmax**(dt/trac_tau)
148  tco_global=zero
149  {^ifoned
150  call mpi_allreduce(tco_mype,tco_global,1,mpi_double_precision,&
151  mpi_max,icomm,ierrmpi)
152  }
153  endif
154  if(.not. associated(phys_trac_after_setdt)) call mpistop("phys_trac_after_setdt not set")
155  ! trac_alfa,tco_global are set only for phys_trac_type=1, should not be a problem when not initialized
156  ! side effect of modifying T_bott from mod_trac -> T_bott sent as param
157  call phys_trac_after_setdt(tco_global,trac_alfa,t_peak, t_bott)
158  end if
159 
160  contains
161 
162  !> compute CFL limited dt (for variable time stepping)
163  subroutine getdt_courant(w,ixI^L,ixO^L,dtnew,dx^D,x,cmax_mype,a2max_mype,cs2max_mype)
167 
168  integer, intent(in) :: ixI^L, ixO^L
169  double precision, intent(in) :: x(ixI^S,1:ndim)
170  double precision, intent(in) :: dx^D
171  double precision, intent(inout) :: w(ixI^S,1:nw), dtnew, cmax_mype, a2max_mype(ndim),cs2max_mype
172 
173  integer :: idims
174  integer :: hxO^L
175  double precision :: courantmax, dxinv(1:ndim), courantmaxtot, courantmaxtots
176  double precision :: cmax(ixI^S), cmaxtot(ixI^S)
177  double precision :: a2max(ndim), cs2max, tco_local, Tmax_local
178 
179  dtnew=bigdouble
180  courantmax=zero
181  courantmaxtot=zero
182  courantmaxtots=zero
183 
184  ! local timestep dt has to be calculated in the
185  ! extended region because of the calculation from the
186  ! div fluxes in mod_finite_volume
187  if(local_timestep) then
188  hxomin^d=ixomin^d-1;
189  hxomax^d=ixomax^d;
190  else
191  hxomin^d=ixomin^d;
192  hxomax^d=ixomax^d;
193  endif
194 
195  if(need_global_a2max) then
196  call phys_get_a2max(w,x,ixi^l,ixo^l,a2max)
197  do idims=1,ndim
198  a2max_mype(idims) = max(a2max_mype(idims),a2max(idims))
199  end do
200  end if
201  if(need_global_cs2max) then
202  call phys_get_cs2max(w,x,ixi^l,ixo^l,cs2max)
203  cs2max_mype = max(cs2max_mype,cs2max)
204  end if
205 
206  if(phys_trac) then
207  call phys_get_tcutoff(ixi^l,ixo^l,w,x,tco_local,tmax_local)
208  {^ifoned tco_mype=max(tco_mype,tco_local) }
209  tmax_mype=max(tmax_mype,tmax_local)
210  end if
211 
212  !if(slab_uniform) then
213  ! ^D&dxinv(^D)=one/dx^D;
214  ! do idims=1,ndim
215  ! call phys_get_cmax(w,x,ixI^L,ixO^L,idims,cmax)
216  ! if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixO^S)))
217  ! if(need_global_a2max) a2max_mype(idims) = max(a2max_mype(idims),a2max(idims))
218  ! cmaxtot(ixO^S)=cmaxtot(ixO^S)+cmax(ixO^S)*dxinv(idims)
219  ! courantmax=max(courantmax,maxval(cmax(ixO^S)*dxinv(idims)))
220  ! courantmaxtot=courantmaxtot+courantmax
221  ! end do
222  !else
223  ! do idims=1,ndim
224  ! call phys_get_cmax(w,x,ixI^L,ixO^L,idims,cmax)
225  ! if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixO^S)))
226  ! if(need_global_a2max) a2max_mype(idims) = max(a2max_mype(idims),a2max(idims))
227  ! tmp(ixO^S)=cmax(ixO^S)/block%ds(ixO^S,idims)
228  ! cmaxtot(ixO^S)=cmaxtot(ixO^S)+tmp(ixO^S)
229  ! courantmax=max(courantmax,maxval(tmp(ixO^S)))
230  ! courantmaxtot=courantmaxtot+courantmax
231  ! end do
232  !end if
233 
234  ! these are also calculated in hxO because of local timestep
235  if(nwaux>0) call phys_get_auxiliary(ixi^l,hxo^l,w,x)
236 
237  select case (type_courant)
238  case (type_maxsum)
239  cmaxtot(hxo^s)=zero
240  if(slab_uniform) then
241  ^d&dxinv(^d)=one/dx^d;
242  do idims=1,ndim
243  call phys_get_cmax(w,x,ixi^l,hxo^l,idims,cmax)
244  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
245  cmaxtot(hxo^s)=cmaxtot(hxo^s)+cmax(hxo^s)*dxinv(idims)
246  end do
247  else
248  do idims=1,ndim
249  call phys_get_cmax(w,x,ixi^l,hxo^l,idims,cmax)
250  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
251  cmaxtot(hxo^s)=cmaxtot(hxo^s)+cmax(hxo^s)/block%ds(hxo^s,idims)
252  end do
253  end if
254  ! courantmaxtots='max(summed c/dx)'
255  courantmaxtots=max(courantmaxtots,maxval(cmaxtot(ixo^s)))
256  if (courantmaxtots>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtots)
257  if(local_timestep) then
258  block%dt(hxo^s) = courantpar/cmaxtot(hxo^s)
259  endif
260 
261  case (type_summax)
262  !TODO this should be mod_input_output?
263  if(local_timestep) then
264  call mpistop("Type courant summax incompatible with local_timestep")
265  endif
266  if(slab_uniform) then
267  ^d&dxinv(^d)=one/dx^d;
268  do idims=1,ndim
269  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
270  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
271  courantmax=max(courantmax,maxval(cmax(ixo^s)*dxinv(idims)))
272  courantmaxtot=courantmaxtot+courantmax
273  end do
274  else
275  do idims=1,ndim
276  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
277  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
278  courantmax=max(courantmax,maxval(cmax(ixo^s)/block%ds(ixo^s,idims)))
279  courantmaxtot=courantmaxtot+courantmax
280  end do
281  end if
282  ! courantmaxtot='summed max(c/dx)'
283  if (courantmaxtot>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtot)
284  case (type_minimum)
285  if(local_timestep) then
286  call mpistop("Type courant not implemented for local_timestep, use maxsum")
287  endif
288  if(slab_uniform) then
289  ^d&dxinv(^d)=one/dx^d;
290  do idims=1,ndim
291  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
292  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
293  courantmax=max(courantmax,maxval(cmax(ixo^s)*dxinv(idims)))
294  end do
295  else
296  do idims=1,ndim
297  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
298  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
299  courantmax=max(courantmax,maxval(cmax(ixo^s)/block%ds(ixo^s,idims)))
300  end do
301  end if
302  ! courantmax='max(c/dx)'
303  if (courantmax>smalldouble) dtnew=min(dtnew,courantpar/courantmax)
304  end select
305  end subroutine getdt_courant
306  end subroutine setdt
307 end module mod_dt
subroutine getdt_courant(w, ixIL, ixOL, dtnew, dxD, x, cmax_mype, a2max_mype, cs2max_mype)
compute CFL limited dt (for variable time stepping)
Definition: mod_dt.t:164
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:11
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.
logical need_global_cs2max
global value for csound speed
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.
double precision cs2max_global
global largest cs2 for hyperbolic thermal conduction
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:62
procedure(sub_get_dt), pointer phys_get_dt
Definition: mod_physics.t:71
procedure(sub_get_tcutoff), pointer phys_get_tcutoff
Definition: mod_physics.t:64
procedure(sub_get_cs2max), pointer phys_get_cs2max
Definition: mod_physics.t:63
procedure(sub_get_auxiliary), pointer phys_get_auxiliary
Definition: mod_physics.t:95
procedure(sub_trac_after_setdt), pointer phys_trac_after_setdt
Definition: mod_physics.t:65
procedure(sub_get_cmax), pointer phys_get_cmax
Definition: mod_physics.t:61
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