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  double precision :: a2max_mype(ndim), tco_mype, tco_global, Tmax_mype, T_bott, T_peak
15  double precision :: trac_alfa, trac_dmax, trac_tau
16 
17  if (dtpar<=zero) then
18  dtmin_mype=bigdouble
19  cmax_mype = zero
20  a2max_mype = zero
21  tco_mype = zero
22  tmax_mype = zero
23  !$OMP PARALLEL DO PRIVATE(igrid,qdtnew,dtnew,dx^D)
24  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
25  dtnew=bigdouble
26  dx^d=rnode(rpdx^d_,igrid);
27  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
28  saveigrid = igrid
29  block=>ps(igrid)
30  block%iw0=0
31 
32  if (nwaux>0) then
33  call phys_get_aux(.true.,ps(igrid)%w,&
34  ps(igrid)%x,ixg^ll,ixm^ll,'setdt')
35  end if
36 
37  call getdt_courant(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,ps(igrid)%x)
38  dtnew=min(dtnew,qdtnew)
39 
40  call phys_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
41  dtnew=min(dtnew,qdtnew)
42 
43  if (associated(usr_get_dt)) then
44  call usr_get_dt(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
45  end if
46 
47  dtnew = min(dtnew,qdtnew)
48  dtmin_mype = min(dtmin_mype,dtnew)
49  dt_grid(igrid) = dtnew
50  end do
51  !$OMP END PARALLEL DO
52  else
53  dtmin_mype=dtpar
54  end if
55 
56  if (dtmin_mype<dtmin) then
57  write(unitterm,*)"Error: Time step too small!", dtmin_mype
58  write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
59  write(unitterm,*)"Lower limit of time step:", dtmin
60  crash=.true.
61  end if
62 
63  if (slowsteps>it-it_init+1) then
64  factor=one-(one-dble(it-it_init+1)/dble(slowsteps))**2
65  dtmin_mype=dtmin_mype*factor
66  end if
67 
68 
69  if(final_dt_reduction)then
70  !if (dtmin_mype>time_max-global_time) then
71  ! write(unitterm,*)"WARNING final timestep artificially reduced!"
72  ! write(unitterm,*)"on processor:", mype, "at time:", global_time," step:", it
73  !endif
74  if(time_max-global_time<=dtmin) then
75  !write(unitterm,*)'Forcing to leave timeloop as time is reached!'
76  final_dt_exit=.true.
77  endif
78  dtmin_mype=min(dtmin_mype,time_max-global_time)
79  endif
80 
81  if (dtpar<=zero) then
82  call mpi_allreduce(dtmin_mype,dt,1,mpi_double_precision,mpi_min, &
83  icomm,ierrmpi)
84  else
85  dt=dtmin_mype
86  end if
87 
88  if(any(dtsave(1:nfile)<bigdouble).or.any(tsave(isavet(1:nfile),1:nfile)<bigdouble))then
89  dtmax = minval(ceiling(global_time/dtsave(1:nfile))*dtsave(1:nfile))-global_time
90  do ifile=1,nfile
91  dtmax = min(tsave(isavet(ifile),ifile)-global_time,dtmax)
92  end do
93  if(dtmax > smalldouble)then
94  dt=min(dt,dtmax)
95  else
96  ! dtmax=0 means dtsave is divisible by global_time
97  dt=min(dt,minval(dtsave(1:nfile)))
98  end if
99  end if
100 
101  if(mype==0) then
102  if(any(dtsave(1:nfile)<dt)) then
103  write(unitterm,*) 'Warning: timesteps: ',dt,' exceeding output intervals ', dtsave(1:nfile)
104  endif
105  endif
106 
107  ! estimate time step of thermal conduction
108  if(associated(phys_getdt_heatconduct)) then
109  dtmin_mype=bigdouble
110  !$OMP PARALLEL DO PRIVATE(igrid,qdtnew,&
111  !$OMP& dx^D) REDUCTION(min:dtmin_mype)
112  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
113  dx^d=rnode(rpdx^d_,igrid);
114  saveigrid = igrid
115  block=>ps(igrid)
116  qdtnew=bigdouble
117  call phys_getdt_heatconduct(ps(igrid)%w,ixg^ll,ixm^ll,qdtnew,dx^d,ps(igrid)%x)
118  dtmin_mype=min(dtmin_mype,qdtnew)
119  end do
120  !$OMP END PARALLEL DO
121  call mpi_allreduce(dtmin_mype,dtnew,1,mpi_double_precision,mpi_min, &
122  icomm,ierrmpi)
123  if(all(flux_scheme=='nul')) dt=min(dt,dtnew)
124  ncycle=ceiling(dt/dtnew)
125  if (ncycle>tc_ncycles) then
126  if(mype==0 .and. .false.) then
127  write(*,*) 'CLF time step is too many times larger than conduction time step',ncycle
128  write(*,*) 'reducing dt to',tc_ncycles,'times of dt_impl!!'
129  endif
130  dt=tc_ncycles*dtnew
131  endif
132  ! get number of sub-steps of supertime stepping (Meyer 2012 MNRAS 422,2102)
133  if(dt/dtnew< 0.5d0) then
134  s=1
135  else if(dt/dtnew< 2.d0) then
136  s=2
137  else
138  s=ceiling((dsqrt(9.d0+8.d0*dt/dtnew)-1.d0)/2.d0)
139  ! only use odd s number
140  s=s/2*2+1
141  endif
142  dt_tc=dt*0.5d0
143  if(mype==0 .and. .false.) write(*,*) 'supertime steps:',s,' normal subcycles:',&
144  ceiling(dt/dtnew/2.d0)
145  endif
146 
147  !$OMP PARALLEL DO PRIVATE(igrid)
148  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
149  dt_grid(igrid)=dt
150  end do
151  !$OMP END PARALLEL DO
152 
153 
154  ! global Lax-Friedrich finite difference flux splitting needs fastest wave-speed
155  ! so does GLM:
156  if(need_global_cmax) call mpi_allreduce(cmax_mype,cmax_global,1,&
157  mpi_double_precision,mpi_max,icomm,ierrmpi)
158  if(need_global_a2max) then
159  call mpi_allreduce(a2max_mype,a2max_global,ndim,&
160  mpi_double_precision,mpi_max,icomm,ierrmpi)
161  end if
162 
163  ! transition region adaptive thermal conduction (Johnston 2019 ApJL, 873, L22)
164  ! transition region adaptive thermal conduction (Johnston 2020 A&A, 635, 168)
165  if(trac) then
166  trac_dmax=0.1d0
167  trac_tau=1.d0/unit_time
168  !> dt or dtnew ?
169  trac_alfa=trac_dmax**(dtnew/trac_tau)
170  {^ifoned
171  call mpi_allreduce(tco_mype,tco_global,1,mpi_double_precision,&
172  mpi_max,icomm,ierrmpi)
173  }
174  call mpi_allreduce(tmax_mype,t_peak,1,mpi_double_precision,&
175  mpi_max,icomm,ierrmpi)
176  ! default lower limit of cutoff temperature
177  t_bott =2.d4/unit_temperature
178  !$OMP PARALLEL DO PRIVATE(igrid)
179  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
180  {^ifoned
181  ps(igrid)%special_values(1)=tco_global
182  }
183  if(ps(igrid)%special_values(1) < t_bott) then
184  ps(igrid)%special_values(1)=t_bott
185  else if(ps(igrid)%special_values(1) > 0.2d0*t_peak) then
186  ps(igrid)%special_values(1)=0.2d0*t_peak
187  end if
188  if(ps(igrid)%special_values(1) .lt. trac_alfa*ps(igrid)%special_values(2)) then
189  ps(igrid)%special_values(1)=trac_alfa*ps(igrid)%special_values(2)
190  end if
191  !> special values(2) to save old tcutoff
192  ps(igrid)%special_values(2)=ps(igrid)%special_values(1)
193  end do
194  !$OMP END PARALLEL DO
195  end if
196 
197  contains
198 
199  !> compute CFL limited dt (for variable time stepping)
200  subroutine getdt_courant(w,ixI^L,ixO^L,dtnew,x)
203 
204  integer, intent(in) :: ixI^L, ixO^L
205  double precision, intent(in) :: x(ixi^s,1:ndim)
206  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
207 
208  integer :: idims
209  double precision :: courantmax, dxinv(1:ndim), courantmaxtot, courantmaxtots
210  double precision :: cmax(ixi^s), cmaxtot(ixi^s), tmp(ixi^s)
211  double precision :: a2max(ndim),tco_local,Tmax_local
212 
213  dtnew=bigdouble
214 
215  courantmax=zero
216  courantmaxtot=zero
217  courantmaxtots=zero
218 
219  ^d&dxinv(^d)=one/dx^d;
220 
221  cmaxtot(ixo^s)=zero
222 
223  if(need_global_a2max) then
224  call phys_get_a2max(w,x,ixi^l,ixo^l,a2max)
225  end if
226  if(trac) then
227  call phys_get_tcutoff(ixi^l,ixo^l,w,x,tco_local,tmax_local)
228  {^ifoned tco_mype=max(tco_mype,tco_local) }
229  tmax_mype=max(tmax_mype,tmax_local)
230  end if
231  do idims=1,ndim
232  call phys_get_cmax(w,x,ixi^l,ixo^l,idims,cmax)
233  if(need_global_cmax) cmax_mype = max(cmax_mype,maxval(cmax(ixo^s)))
234  if(need_global_a2max) a2max_mype = max(a2max_mype,a2max(idims))
235  if(slab_uniform) then
236  cmaxtot(ixo^s)=cmaxtot(ixo^s)+cmax(ixo^s)*dxinv(idims)
237  courantmax=max(courantmax,maxval(cmax(ixo^s)*dxinv(idims)))
238  else
239  tmp(ixo^s)=cmax(ixo^s)/block%ds(ixo^s,idims)
240  cmaxtot(ixo^s)=cmaxtot(ixo^s)+tmp(ixo^s)
241  courantmax=max(courantmax,maxval(tmp(ixo^s)))
242  end if
243  courantmaxtot=courantmaxtot+courantmax
244  end do
245 
246  select case (typecourant)
247  case ('minimum')
248  ! courantmax='max(c/dx)'
249  if (courantmax>smalldouble) dtnew=min(dtnew,courantpar/courantmax)
250  case ('summax')
251  ! courantmaxtot='summed max(c/dx)'
252  if (courantmaxtot>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtot)
253  case ('maxsum')
254  ! courantmaxtots='max(summed c/dx)'
255  courantmaxtots=max(courantmaxtots,maxval(cmaxtot(ixo^s)))
256  if (courantmaxtots>smalldouble) dtnew=min(dtnew,courantpar/courantmaxtots)
257  case default
258  write(unitterm,*)'Unknown typecourant=',typecourant
259  call mpistop("Error from getdt_courant: no such typecourant!")
260  end select
261 
262  end subroutine getdt_courant
263 
264 end subroutine setdt
This module contains definitions of global parameters and variables and some generic functions/subrou...
logical need_global_a2max
global value for schmid scheme
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, dimension(ndim) a2max_global
global largest a2 for schmid scheme
procedure(sub_get_a2max), pointer phys_get_a2max
Definition: mod_physics.t:66
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(sub_get_tcutoff), pointer phys_get_tcutoff
Definition: mod_physics.t:67
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
logical final_dt_reduction
If true, allow final dt reduction for matching time_max on output.
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:72
double precision, dimension(nfile) dtsave
Repeatedly save output of type N when dtsave(N) simulation time has passed.
double precision unit_temperature
Physical scaling factor for temperature.
subroutine getdt_courant(w, ixIL, ixOL, dtnew, x)
compute CFL limited dt (for variable time stepping)
Definition: setdt.t:201
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:65
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:77
double precision unit_time
Physical scaling factor for time.
character(len=std_len) typecourant
How to compute the CFL-limited time step.
logical trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
integer it
Number of time steps taken.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
logical final_dt_exit
Force timeloop exit when final dt < dtmin.