MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_supertimestepping.t
Go to the documentation of this file.
1 !> Generic supertimestepping method
2 !> 1) in amrvac.par in sts_list set the following parameters which have the default values:
3 !> sts_dtpar=0.9,sts_ncycles=1000,sts_method=1,sourcetype_sts=2
4 !> These parametes are general for all the methods added TODO: check if there is any need
5 !> to have terms implemented with different sets of parameters, and these cannot be general anymore
6 !> 2) then add programatically in the code a term with the subroutine
7 !> add_sts_method
8 !> This method takes as parameters a function which calculated the explicit timestep
9 !> associated with the term, a subroutine which sets the source term
10 !> types for the BC and the BC are generated from the variables startVar:endVar
11 !> flux conservation (fixconserve) is done for the variables specified by ixChangeStart, ixChangeN, ixChangeFixC
12 !> The following two steps are done in this way as in fortran it is not allowed to pass null function pointers as parameters:
13 !> 3)in order to to have hooks before_first_cycle, after_last_cycle (e.g. conversion from e_tot to e_int before first sts cycle
14 !> and back from e_int to e_tot after the last STS cycle for the thermal conductivity module) add them just afterwards with the subroutine
15 !> set_conversion_methods_to_head
16 !> 4) to add the hook for error handling (e.g check small values in the thermal conductivity module )
17 !> call set_error_handling_to_head which takes as parameter a subroutine
18 !> the error handling subroutine is called before setting BC
20  use mod_geometry
21  implicit none
22  private
23 
24  public :: is_sts_initialized
25  public :: sts_init
27  public :: sts_add_source
28  public :: set_dt_sts_ncycles
30 
31  ! input parameters from parameter file
32  !> the coefficient that multiplies the sts dt
33  double precision :: sts_dtpar=0.8d0
34  !> the maximum number of subcycles
35  integer :: sts_ncycles=1000
36  integer :: sts_method = 1
37  integer, parameter :: sourcetype_sts_prior =0
38  integer, parameter :: sourcetype_sts_after =1
39  integer, parameter :: sourcetype_sts_split =2
41 
42  !The following is used only for method 2, not input parameter TODO check if we want as input parameter
43  double precision,parameter :: nu_sts = 0.5d0
44  !> Whether to conserve fluxes at the current partial step
45  logical :: fix_conserve_at_step = .true.
46  logical :: sts_initialized = .false.
47 
48  abstract interface
49 
50  !>interface for setting sources in the derived type
51  subroutine subr1(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
53  integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
54  double precision, intent(in) :: x(ixi^s,1:ndim)
55  double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
56  double precision, intent(in) :: my_dt
57  logical, intent(in) :: fix_conserve_at_step
58  end subroutine subr1
59 
60  !>interface for the function which gets the timestep in dtnew in the derived type
61  function subr2(w,ixG^L,ix^L,dx^D,x) result(dtnew)
63  integer, intent(in) :: ixG^L, ix^L
64  double precision, intent(in) :: dx^D, x(ixG^S,1:ndim)
65  double precision, intent(in) :: w(ixG^S,1:nw)
66  double precision :: dtnew
67  end function subr2
68 
69  !>interface for error handling subroutine in the derived type
70  subroutine subr_e(w, x, ixI^L, ixO^L, step)
73  integer, intent(in) :: ixI^L,ixO^L
74  double precision, intent(inout) :: w(ixI^S,1:nw)
75  double precision, intent(in) :: x(ixI^S,1:ndim)
76  integer, intent(in) :: step
77  end subroutine subr_e
78 
79  !>interface for the subroutines before_first_cycle and after_last_cycle in the derived type
80  subroutine subr5(ixI^L, ixO^L, w, x)
82  integer, intent(in) :: ixI^L, ixO^L
83  double precision, intent(in) :: x(ixI^S,1:ndim)
84  double precision, intent(inout) :: w(ixI^S,1:nw)
85  end subroutine subr5
86 
87  !>for the subroutines in this module, which do not depend on the term, but
88  !>on the parameter sts_method = 1/2 in the parameter file
89  !>sts_add_source
90  subroutine subr3(dt)
91  double precision,intent(in) :: dt
92  end subroutine subr3
93 
94  !>sts_get_ncycles
95  function subr4(dt,dtnew,dt_modified) result(s)
96  double precision,intent(in) :: dtnew
97  double precision,intent(inout) :: dt
98  logical,intent(inout) :: dt_modified
99  integer :: s
100  end function subr4
101 
102  end interface
103 
104  type sts_term
105 
106  procedure(subr1), pointer, nopass :: sts_set_sources
107  procedure(subr2), pointer, nopass :: sts_getdt
108  procedure(subr5), pointer, nopass :: sts_before_first_cycle, sts_after_last_cycle
109  procedure(subr_e), pointer, nopass :: sts_handle_errors
110  type(sts_term), pointer :: next
111  double precision :: dt_expl
112  integer, public :: s
113  logical :: types_initialized
114  logical :: evolve_magnetic_field
115 
116  !>types used for send/recv ghosts, see mod_ghostcells_update
117  integer, dimension(-1:2^D&,-1:1^D&) :: type_send_srl_sts_1, type_recv_srl_sts_1
118  integer, dimension(-1:1^D&,-1:1^D&) :: type_send_r_sts_1
119  integer, dimension(-1:1^D&, 0:3^D&) :: type_recv_r_sts_1
120  integer, dimension(-1:1^D&, 0:3^D&) :: type_recv_p_sts_1, type_send_p_sts_1
121 
122  integer, dimension(-1:2^D&,-1:1^D&) :: type_send_srl_sts_2, type_recv_srl_sts_2
123  integer, dimension(-1:1^D&,-1:1^D&) :: type_send_r_sts_2
124  integer, dimension(-1:1^D&, 0:3^D&) :: type_recv_r_sts_2
125  integer, dimension(-1:1^D&, 0:3^D&) :: type_recv_p_sts_2, type_send_p_sts_2
126 
127  integer :: startVar
128  integer :: endVar
129  !> number of flux involved in STS update
130  integer :: nflux
131  integer :: startwbc
132  integer :: nwbc
133 
134  end type sts_term
135 
136  type(sts_term), pointer :: head_sts_terms
137  !The following two subroutine/function pointers
138  !make the difference between the two STS methods implemented
139  procedure(subr3), pointer :: sts_add_source
140  procedure(subr4), pointer :: sts_get_ncycles
141 
142 contains
143 
144  !> Initialize sts module
145  subroutine sts_init()
147  use mod_physics
148  if(.not. sts_initialized) then
149  nullify(head_sts_terms)
150  call sts_params_read(par_files)
151  sts_dtpar=sts_dtpar/dble(ndim)
152  sts_initialized = .true.
153  if(sts_method .eq. 1) then
155  sts_get_ncycles => sts_get_ncycles1
156  else if(sts_method .eq. 2) then
157  sts_add_source => sts_add_source2
158  sts_get_ncycles => sts_get_ncycles2
159  else
160  call mpistop("Unknown sts method")
161  end if
162  endif
163 
164  end subroutine sts_init
165 
166  pure function is_sts_initialized() result(res)
167  logical :: res
168  if (sts_initialized) then
169  res = associated(head_sts_terms)
170  else
171  res = .false.
172  endif
173  end function is_sts_initialized
174 
175  !> Read module parameters from par file
176  subroutine sts_params_read(files)
177  use mod_global_parameters, only: unitpar
178  character(len=*), intent(in) :: files(:)
179  integer :: n
180 
181  namelist /sts_list/ sts_dtpar,sts_ncycles,sts_method,sourcetype_sts
182 
183  do n = 1, size(files)
184  open(unitpar, file=trim(files(n)), status="old")
185  read(unitpar, sts_list, end=111)
186 111 close(unitpar)
187  end do
188 
189  end subroutine sts_params_read
190 
191  !> subroutine which added programatically a term to be calculated using STS
192  !> Params:
193  !> sts_getdt function calculates the explicit timestep for this term
194  !> sts_set_sources subroutine sets the source term
195  !> startVar, endVar, nflux indices of start, end, and number of the variables that need fix conservation
196  !> startwbc, nwbc indices of start and number of the variables that need ghost cell exchange
197  !> These terms implemented by an element of the derived type sts_term are put in a linked list
198  subroutine add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
201 
202  integer, intent(in) :: startvar, nflux, startwbc, nwbc
203  logical, intent(in) :: evolve_b
204 
205  interface
206 
207  subroutine sts_set_sources(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
209  use mod_fix_conserve
210  integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
211  double precision, intent(in) :: x(ixi^s,1:ndim)
212  double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
213  double precision, intent(in) :: my_dt
214  logical, intent(in) :: fix_conserve_at_step
215  end subroutine sts_set_sources
216 
217  function sts_getdt(w,ixG^L,ix^L,dx^D,x) result(dtnew)
219  integer, intent(in) :: ixg^l, ix^l
220  double precision, intent(in) :: dx^d, x(ixg^s,1:ndim)
221  double precision, intent(in) :: w(ixg^s,1:nw)
222  double precision :: dtnew
223  end function sts_getdt
224 
225  end interface
226 
227  type(sts_term), pointer :: temp
228  allocate(temp)
229 
230  temp%sts_getdt => sts_getdt
231  temp%sts_set_sources => sts_set_sources
232  temp%sts_before_first_cycle => null()
233  temp%sts_after_last_cycle => null()
234  temp%sts_handle_errors => null()
235  temp%startVar = startvar
236  temp%endVar= startvar+nflux-1
237  temp%nflux = nflux
238  temp%startwbc = startwbc
239  temp%nwbc = nwbc
240  temp%types_initialized = .false.
241  temp%evolve_magnetic_field=evolve_b
242 
243  temp%next => head_sts_terms
244  head_sts_terms => temp
245 
246  end subroutine add_sts_method
247 
248  !> Set the hooks called before the first cycle and after the last cycle in the STS update
249  !> This method should be called after add_sts_method. The hooks are added to the last term added with this subroutine
250  !> Params: sts_before_first_cycle, sts_after_last_cycle subroutines which implement the hooks called before first cycle and after last cycle
251  subroutine set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
252  interface
253  subroutine sts_before_first_cycle(ixI^L, ixO^L, w, x)
255  integer, intent(in) :: ixi^l, ixo^l
256  double precision, intent(in) :: x(ixi^s,1:ndim)
257  double precision, intent(inout) :: w(ixi^s,1:nw)
258  end subroutine sts_before_first_cycle
259 
260  subroutine sts_after_last_cycle(ixI^L, ixO^L, w, x)
262  integer, intent(in) :: ixi^l, ixo^l
263  double precision, intent(in) :: x(ixi^s,1:ndim)
264  double precision, intent(inout) :: w(ixi^s,1:nw)
265  end subroutine sts_after_last_cycle
266  end interface
267 
268  head_sts_terms%sts_before_first_cycle => sts_before_first_cycle
269  head_sts_terms%sts_after_last_cycle => sts_after_last_cycle
270 
271  end subroutine set_conversion_methods_to_head
272 
273  !> Set the hook of error handling in the STS update. This method is called before updating the BC.
274  !> This method should be called after add_sts_method. The hook is added to the last term added with this subroutine.
275  !> Param: sts_error_handing the subroutine which handles the errors
276  subroutine set_error_handling_to_head(sts_error_handling)
277  interface
278  subroutine sts_error_handling(w, x, ixI^L, ixO^L, step)
280  use mod_small_values
281  integer, intent(in) :: ixi^l,ixo^l
282  double precision, intent(inout) :: w(ixi^s,1:nw)
283  double precision, intent(in) :: x(ixi^s,1:ndim)
284  integer, intent(in) :: step
285  end subroutine sts_error_handling
286  end interface
287  head_sts_terms%sts_handle_errors => sts_error_handling
288 
289  end subroutine set_error_handling_to_head
290 
291  !> method used to set the number of cycles for the STS1 method
292  function sts_get_ncycles1(dt,dtnew,dt_modified) result(is)
293  double precision,intent(in) :: dtnew
294  double precision,intent(inout) :: dt
295  logical,intent(inout) :: dt_modified
296  integer :: is
297 
298  double precision :: ss
299 
300  !!ss is now limit of dt because of sts_ncycles
301  ss = dtnew*((2.d0*sts_ncycles+1)**2-9.d0)/16.d0
302  if(dt>ss) then
303  dt_modified = .true.
304  dt = ss
305  is = sts_ncycles
306  else
307  ss = dt/dtnew
308  ! get number of sub-steps of supertime stepping (Meyer 2012 MNRAS 422,2102)
309  if(ss .le. 1.d0) then
310  is=1
311  else
312  is=ceiling((dsqrt(9.d0+16.d0*ss)-1.d0)*0.5d0)
313  is=is/2*2+1
314  end if
315  end if
316 
317  end function sts_get_ncycles1
318 
319  !> method used to set the number of cycles for the STS2 method
320  function sts_get_ncycles2(dt,dtnew,dt_modified) result(is)
321  double precision,intent(in) :: dtnew
322  double precision,intent(inout) :: dt
323  logical,intent(inout) :: dt_modified
324  integer :: is
325 
326  double precision :: ss,rr
327  integer:: ncycles
328 
329  rr = dt/dtnew
330  !print*, dt, " --DTEXPL-- ", dtnew, ", rr ",rr
331  ncycles = sts_ncycles
332  !print*, "NCYCLES BEFORE ",ncycles
333  ss=sum_chev(nu_sts,ncycles,rr)
334  !print*, "NCYCLES AFTER ",ncycles
335  is = ncycles
336  !print*, "SUMCHEV ", ss, " NCYCLES ", is
337  if(ss < rr) then
338  dt_modified = .true.
339  dt = ss * dtnew
340  endif
341 
342  end function sts_get_ncycles2
343 
344  !> This sets the explicit dt and calculates the number of cycles for each of the terms implemented with STS.
345  function set_dt_sts_ncycles(my_dt) result(dt_modified)
347 
348  double precision,intent(inout) :: my_dt
349  double precision :: my_dt1
350  logical :: dt_modified, dt_modified1, dt_modified2
351 
352  integer:: iigrid, igrid, ncycles
353  double precision :: dtnew,dtmin_mype
354  double precision :: dx^d, ss
355  type(sts_term), pointer :: temp,oldtemp
356  nullify(oldtemp)
357  temp => head_sts_terms
358  dt_modified = .false.
359  do while(associated(temp))
360  dt_modified2 = .false.
361  dtmin_mype=bigdouble
362  !$OMP PARALLEL DO PRIVATE(igrid,dx^D) REDUCTION(min:dtmin_mype)
363  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
364  ! maybe the following global variables are needed in get_dt!
365  ! next few lines ensure correct usage of routines like divvector etc
366  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
367  block=>ps(igrid)
368  ! end maybe the following global variables are needed in get_dt!!!!!!!
369  dx^d=rnode(rpdx^d_,igrid);
370  dtmin_mype=min(dtmin_mype, sts_dtpar * temp%sts_getdt(ps(igrid)%w,ixg^ll,ixm^ll,dx^d,ps(igrid)%x))
371  end do
372  !$OMP END PARALLEL DO
373  call mpi_allreduce(dtmin_mype,dtnew,1,mpi_double_precision,mpi_min,icomm,ierrmpi)
374  temp%s = sts_get_ncycles(my_dt,dtnew,dt_modified2)
375 
376  !print*, "NCYCLES ", temp%s, dt_modified2, my_dt, dtnew
377  temp%dt_expl = dtnew
378 
379  ! Note that as for some term it may happen that the dt is modified: it may be reduced if the
380  ! number of cycles is overpassed, the list has to be reiterated to update ncycles for previous
381  ! terms which did not modify dt TODO add pointer to previous and loop backward to update
382  if(dt_modified2) then
383  dt_modified = .true.
384  !reiterate all the other sts elements and recalculate s
385  oldtemp => head_sts_terms
386  my_dt1 = my_dt
387  dt_modified1 = .false.
388  do while(.not. associated(oldtemp,temp))
389  oldtemp%s = sts_get_ncycles(my_dt1,oldtemp%dt_expl,dt_modified1)
390  !check dt is not modified again, and this should not happen, except for bug in sts_get_ncycles1,2
391  if(dt_modified1) call mpistop("sts dt modified twice")
392  oldtemp=>oldtemp%next
393  end do
394  end if
395  temp=>temp%next
396 
397  end do
398 
399  end function set_dt_sts_ncycles
400 
401  pure FUNCTION chev(j,nu,N)
402  use mod_constants
403 
404  double precision, INTENT(IN) :: nu
405  INTEGER, INTENT(IN) :: j, n
406  double precision :: chev
407 
408  chev = 1d0 / ((-1d0 + nu)*cos(((2d0*j - 1d0) / n)* (dpi/2d0)) + 1d0 + nu)
409 
410  END FUNCTION chev
411 
412  FUNCTION sum_chev(nu,N,limMax)
413  double precision, intent(in) :: nu,limmax
414  integer, intent(inout) :: n
415  double precision :: sum_chev, tmp
416 
417  integer :: j
418 
419  j=1
420  sum_chev = 0d0
421  do while (j < n .and. sum_chev < limmax)
422  sum_chev = sum_chev + chev(j,nu,n)
423  j=j+1
424  enddo
425  n=j-1
426  END FUNCTION sum_chev
427 
428  !TODO the following not used
429 ! PURE FUNCTION total_chev(nu,N)
430 ! double precision, INTENT(IN) :: nu
431 ! INTEGER, INTENT(IN) :: N
432 ! double precision :: total_chev
433 !
434 ! total_chev = N/(2d0*dsqrt(nu)) * ( (1d0 + dsqrt(nu))**(2d0*N) - (1d0 - dsqrt(nu))**(2d0*N) ) / &
435 ! ( (1d0 + dsqrt(nu))**(2d0*N) + (1d0 - dsqrt(nu))**(2d0*N) )
436 !
437 ! END FUNCTION total_chev
438 
439  !> Iterates all the terms implemented with STS and adds the sources
440  !> STS method 2 implementation
441  subroutine sts_add_source2(my_dt)
442  ! Turlough Downes 2006,2007
445  use mod_fix_conserve
446  use mod_physics
447 
448  double precision, intent(in) :: my_dt
449  double precision, allocatable :: bj(:)
450  double precision :: sumbj,dtj
451 
452  integer:: iigrid, igrid, j, ixc^l
453  logical :: stagger_flag=.false., prolong_flag=.false., coarsen_flag=.false.
454  type(sts_term), pointer :: temp
455 
456  ! do not fill physical boundary conditions
457  bcphys=.false.
458 
459  fix_conserve_at_step = time_advance .and. levmax>levmin
460 
461  temp => head_sts_terms
462  do while(associated(temp))
463 
464  if(.not.temp%evolve_magnetic_field) then
465  ! not do fix conserve and getbc for staggered values
466  stagger_flag=stagger_grid
467  stagger_grid=.false.
468  else if(stagger_grid) then
469  ixcmax^d=ixmhi^d;
470  ixcmin^d=ixmlo^d-1;
471  end if
472 
473  call init_comm_fix_conserve(1,ndim,temp%nflux)
474 
475  if(associated(temp%sts_before_first_cycle)) then
476  prolong_flag=prolongprimitive
477  coarsen_flag=coarsenprimitive
478  prolongprimitive=.false.
479  coarsenprimitive=.false.
480  do iigrid=1,igridstail; igrid=igrids(iigrid);
481  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
482  block=>ps(igrid)
483  call temp%sts_before_first_cycle(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
484  end do
485  end if
486 
487  allocate(bj(1:temp%s))
488  do j=1,temp%s
489  bj(j) = chev(j,nu_sts,sts_ncycles)
490  end do
491 
492  type_send_srl=>temp%type_send_srl_sts_1
493  type_recv_srl=>temp%type_recv_srl_sts_1
494  type_send_r=>temp%type_send_r_sts_1
495  type_recv_r=>temp%type_recv_r_sts_1
496  type_send_p=>temp%type_send_p_sts_1
497  type_recv_p=>temp%type_recv_p_sts_1
498 
499  if(.not. temp%types_initialized) then
500  call create_bc_mpi_datatype(temp%startwbc,temp%nwbc)
501  if(temp%nflux>temp%nwbc) then
502  ! prepare types for the changed no-need-ghost-update variables in the last getbc
503  type_send_srl=>temp%type_send_srl_sts_2
504  type_recv_srl=>temp%type_recv_srl_sts_2
505  type_send_r=>temp%type_send_r_sts_2
506  type_recv_r=>temp%type_recv_r_sts_2
507  type_send_p=>temp%type_send_p_sts_2
508  type_recv_p=>temp%type_recv_p_sts_2
509  call create_bc_mpi_datatype(temp%startVar,temp%nflux)
510  type_send_srl=>temp%type_send_srl_sts_1
511  type_recv_srl=>temp%type_recv_srl_sts_1
512  type_send_r=>temp%type_send_r_sts_1
513  type_recv_r=>temp%type_recv_r_sts_1
514  type_send_p=>temp%type_send_p_sts_1
515  type_recv_p=>temp%type_recv_p_sts_1
516  end if
517  temp%types_initialized = .true.
518  end if
519 
520  sumbj=0.d0
521  do j=1,temp%s
522  if(j .eq. temp%s .and. (sumbj + bj(j)) * temp%dt_expl > my_dt) then
523  dtj = my_dt - sumbj * temp%dt_expl
524  else
525  dtj = bj(j)* temp%dt_expl
526  end if
527  sumbj = sumbj + bj(j)
528  if(stagger_grid) then
529  !$OMP PARALLEL DO PRIVATE(igrid)
530  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
531  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
532  block=>ps(igrid)
533  call temp%sts_set_sources(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,ps1(igrid)%w,fix_conserve_at_step,dtj,igrid,temp%nflux)
534  if(temp%nflux>ndir) then
535  ps(igrid)%w(ixm^t,temp%startVar)=ps(igrid)%w(ixm^t,temp%startVar)+dtj*ps1(igrid)%w(ixm^t,temp%startVar)
536  end if
537  ps(igrid)%ws(ixc^s,1:nws)=ps(igrid)%ws(ixc^s,1:nws)+dtj*ps1(igrid)%w(ixc^s,iw_mag(1:nws))
538  call phys_face_to_center(ixm^ll,ps(igrid))
539  end do
540  !$OMP END PARALLEL DO
541  else
542  !$OMP PARALLEL DO PRIVATE(igrid)
543  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
544  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
545  block=>ps(igrid)
546  call temp%sts_set_sources(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,ps1(igrid)%w,fix_conserve_at_step,dtj,igrid,temp%nflux)
547  ps(igrid)%w(ixm^t,temp%startVar:temp%endVar)=ps(igrid)%w(ixm^t,temp%startVar:temp%endVar)+&
548  dtj*ps1(igrid)%w(ixm^t,temp%startVar:temp%endVar)
549  end do
550  !$OMP END PARALLEL DO
551  end if
552  !fix conserve the fluxes set in the STS method
553  if(fix_conserve_at_step) then
554  call recvflux(1,ndim)
555  call sendflux(1,ndim)
556  call fix_conserve(ps,1,ndim,temp%startVar,temp%nflux)
557  if(stagger_grid) then
558  call fix_edges(ps,1,ndim)
559  ! fill the cell-center values from the updated staggered variables
560  !$OMP PARALLEL DO PRIVATE(igrid)
561  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
562  call phys_face_to_center(ixg^ll,ps(igrid))
563  end do
564  !$OMP END PARALLEL DO
565  end if
566  end if
567  if(associated(temp%sts_handle_errors)) then
568  !$OMP PARALLEL DO PRIVATE(igrid)
569  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
570  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
571  block=>ps(igrid)
572  call temp%sts_handle_errors(ps(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,j)
573  end do
574  !$OMP END PARALLEL DO
575  end if
576 
577  if(temp%nflux>temp%nwbc.and.temp%s==j) then
578  ! include the changed no-need-ghost-update variables in the last getbc
579  type_send_srl=>temp%type_send_srl_sts_2
580  type_recv_srl=>temp%type_recv_srl_sts_2
581  type_send_r=>temp%type_send_r_sts_2
582  type_recv_r=>temp%type_recv_r_sts_2
583  type_send_p=>temp%type_send_p_sts_2
584  type_recv_p=>temp%type_recv_p_sts_2
585  call getbc(global_time,0.d0,ps,temp%startVar,temp%nflux)
586  else
587  call getbc(global_time,0.d0,ps,temp%startwbc,temp%nwbc)
588  end if
589  end do
590 
591  if(associated(temp%sts_after_last_cycle)) then
592  do iigrid=1,igridstail; igrid=igrids(iigrid);
593  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
594  block=>ps(igrid)
595  call temp%sts_after_last_cycle(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
596  end do
597  prolongprimitive = prolong_flag
598  coarsenprimitive = coarsen_flag
599  end if
600  deallocate(bj)
601 
602  if(.not.temp%evolve_magnetic_field) then
603  ! restore stagger_grid value
604  stagger_grid=stagger_flag
605  end if
606 
607  temp=>temp%next
608  end do
609 
610  if(associated(head_sts_terms)) then
611  ! point bc mpi data type back to full type for (M)HD
618  end if
619 
620  bcphys=.true.
621 
622  if(phys_partial_ionization) then
623  ! update temperature variable in w
624  !$OMP PARALLEL DO PRIVATE(igrid)
625  do iigrid=1,igridstail; igrid=igrids(iigrid);
626  call phys_update_temperature(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%w,ps(igrid)%x)
627  end do
628  !$OMP END PARALLEL DO
629  end if
630 
631  end subroutine sts_add_source2
632 
633  !> Iterates all the terms implemented with STS and adds the sources
634  !> STS method 1 implementation
635  subroutine sts_add_source1(my_dt)
636  ! Meyer 2012 MNRAS 422,2102
639  use mod_fix_conserve
640  use mod_physics
642 
643  double precision, intent(in) :: my_dt
644  double precision :: dtj
645  double precision :: omega1,cmu,cmut,cnu,cnut,one_mu_nu
646  double precision, allocatable :: bj(:)
647  integer:: iigrid, igrid, j, ixC^L, ixGext^L
648  logical :: evenstep, stagger_flag=.false., prolong_flag=.false., coarsen_flag=.false., total_energy_flag=.true.
649  type(sts_term), pointer :: temp
650  type(state), dimension(:), pointer :: tmpPs1, tmpPs2
651 
652  ! do not fill physical boundary conditions
653  bcphys=.false.
654 
655  fix_conserve_at_step = time_advance .and. levmax>levmin
656 
657  temp => head_sts_terms
658  do while(associated(temp))
659 
660  if(.not.temp%evolve_magnetic_field) then
661  ! not do fix conserve and getbc for staggered values
662  stagger_flag=stagger_grid
663  stagger_grid=.false.
664  else if(stagger_grid) then
665  ixcmax^d=ixmhi^d;
666  ixcmin^d=ixmlo^d-1;
667  end if
668 
669  call init_comm_fix_conserve(1,ndim,temp%nflux)
670 
671  if(associated(temp%sts_before_first_cycle)) then
672  prolong_flag = prolongprimitive
673  coarsen_flag = coarsenprimitive
674  prolongprimitive = .false.
675  coarsenprimitive = .false.
676  total_energy_flag=phys_total_energy
677  phys_total_energy=.false.
678  !$OMP PARALLEL DO PRIVATE(igrid)
679  do iigrid=1,igridstail; igrid=igrids(iigrid);
680  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
681  block=>ps(igrid)
682  call temp%sts_before_first_cycle(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
683  if(.not. allocated(ps2(igrid)%w)) allocate(ps2(igrid)%w(ixg^t,1:nw))
684  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
685  if(.not. allocated(ps4(igrid)%w)) allocate(ps4(igrid)%w(ixg^t,1:nw))
686  ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
687  ps2(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
688  end do
689  !$OMP END PARALLEL DO
690  else
691  if(stagger_grid) then
692  ixgext^l=ixg^ll^ladd1;
693  !$OMP PARALLEL DO PRIVATE(igrid)
694  do iigrid=1,igridstail; igrid=igrids(iigrid);
695  if(.not. allocated(ps2(igrid)%w)) then
696  call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
697  end if
698  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
699  if(.not. allocated(ps4(igrid)%w)) allocate(ps4(igrid)%w(ixg^t,1:nw))
700  ps1(igrid)%w=ps(igrid)%w
701  ps2(igrid)%w=ps(igrid)%w
702  ps1(igrid)%ws=ps(igrid)%ws
703  ps2(igrid)%ws=ps(igrid)%ws
704  end do
705  !$OMP END PARALLEL DO
706  else
707  !$OMP PARALLEL DO PRIVATE(igrid)
708  do iigrid=1,igridstail; igrid=igrids(iigrid);
709  if(.not. allocated(ps2(igrid)%w)) allocate(ps2(igrid)%w(ixg^t,1:nw))
710  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
711  if(.not. allocated(ps4(igrid)%w)) allocate(ps4(igrid)%w(ixg^t,1:nw))
712  ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
713  ps2(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
714  end do
715  !$OMP END PARALLEL DO
716  end if
717  end if
718 
719  allocate(bj(0:temp%s))
720  bj(0)=1.d0/3.d0
721  bj(1)=bj(0)
722  if(temp%s>1) then
723  omega1=4.d0/dble(temp%s**2+temp%s-2)
724  cmut=omega1/3.d0
725  else
726  omega1=0.d0
727  cmut=1.d0
728  end if
729 
730  type_send_srl=>temp%type_send_srl_sts_1
731  type_recv_srl=>temp%type_recv_srl_sts_1
732  type_send_r=>temp%type_send_r_sts_1
733  type_recv_r=>temp%type_recv_r_sts_1
734  type_send_p=>temp%type_send_p_sts_1
735  type_recv_p=>temp%type_recv_p_sts_1
736 
737  if(.not. temp%types_initialized) then
738  call create_bc_mpi_datatype(temp%startwbc,temp%nwbc)
739  if(temp%nflux>temp%nwbc) then
740  ! prepare types for the changed no-need-ghost-update variables in the last getbc
741  type_send_srl=>temp%type_send_srl_sts_2
742  type_recv_srl=>temp%type_recv_srl_sts_2
743  type_send_r=>temp%type_send_r_sts_2
744  type_recv_r=>temp%type_recv_r_sts_2
745  type_send_p=>temp%type_send_p_sts_2
746  type_recv_p=>temp%type_recv_p_sts_2
747  call create_bc_mpi_datatype(temp%startVar,temp%nflux)
748  type_send_srl=>temp%type_send_srl_sts_1
749  type_recv_srl=>temp%type_recv_srl_sts_1
750  type_send_r=>temp%type_send_r_sts_1
751  type_recv_r=>temp%type_recv_r_sts_1
752  type_send_p=>temp%type_send_p_sts_1
753  type_recv_p=>temp%type_recv_p_sts_1
754  end if
755  temp%types_initialized = .true.
756  end if
757  dtj = cmut*my_dt
758  if(stagger_grid) then
759  !$OMP PARALLEL DO PRIVATE(igrid)
760  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
761  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
762  block=>ps(igrid)
763  ps4(igrid)%w=zero
764  call temp%sts_set_sources(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,ps4(igrid)%w,fix_conserve_at_step,dtj,igrid,temp%nflux)
765  !!!eq solved: dU/dt = S, ps3 is stored S^n
766  ps3(igrid)%w(ixc^s,temp%startVar:temp%endVar) = my_dt * ps4(igrid)%w(ixc^s,temp%startVar:temp%endVar)
767  if(temp%nflux>ndir) then
768  ps1(igrid)%w(ixm^t,temp%startVar) = ps1(igrid)%w(ixm^t,temp%startVar) + cmut * ps3(igrid)%w(ixm^t,temp%startVar)
769  end if
770  ps1(igrid)%ws(ixc^s,1:nws) = ps1(igrid)%ws(ixc^s,1:nws) + cmut * ps3(igrid)%w(ixc^s,iw_mag(1:nws))
771  call phys_face_to_center(ixm^ll,ps1(igrid))
772  end do
773  !$OMP END PARALLEL DO
774  else
775  !$OMP PARALLEL DO PRIVATE(igrid)
776  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
777  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
778  block=>ps(igrid)
779  call temp%sts_set_sources(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,ps4(igrid)%w,fix_conserve_at_step,dtj,igrid,temp%nflux)
780  !!!eq solved: dU/dt = S, ps3 is stored S^n
781  ps3(igrid)%w(ixm^t,temp%startVar:temp%endVar) = my_dt * ps4(igrid)%w(ixm^t,temp%startVar:temp%endVar)
782  ps1(igrid)%w(ixm^t,temp%startVar:temp%endVar) = ps1(igrid)%w(ixm^t,temp%startVar:temp%endVar) + &
783  cmut * ps3(igrid)%w(ixm^t,temp%startVar:temp%endVar)
784  end do
785  !$OMP END PARALLEL DO
786  end if
787  if(fix_conserve_at_step) then
788  call recvflux(1,ndim)
789  call sendflux(1,ndim)
790  call fix_conserve(ps1,1,ndim,temp%startVar,temp%nflux)
791  if(stagger_grid) then
792  call fix_edges(ps1,1,ndim)
793  ! fill the cell-center values from the updated staggered variables
794  !$OMP PARALLEL DO PRIVATE(igrid)
795  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
796  call phys_face_to_center(ixg^ll,ps1(igrid))
797  end do
798  !$OMP END PARALLEL DO
799  end if
800  end if
801  ! fix conservation of AMR grid by replacing flux from finer neighbors
802  if(associated(temp%sts_handle_errors)) then
803  !$OMP PARALLEL DO PRIVATE(igrid)
804  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
805  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
806  block=>ps(igrid)
807  call temp%sts_handle_errors(ps1(igrid)%w,ps1(igrid)%x,ixg^ll,ixm^ll,1)
808  end do
809  !$OMP END PARALLEL DO
810  end if
811  if(temp%nflux>temp%nwbc.and.temp%s==1) then
812  ! include the changed no-need-ghost-update variables in the last getbc
813  type_send_srl=>temp%type_send_srl_sts_2
814  type_recv_srl=>temp%type_recv_srl_sts_2
815  type_send_r=>temp%type_send_r_sts_2
816  type_recv_r=>temp%type_recv_r_sts_2
817  type_send_p=>temp%type_send_p_sts_2
818  type_recv_p=>temp%type_recv_p_sts_2
819  call getbc(global_time,0.d0,ps1,temp%startVar,temp%nflux)
820  else
821  call getbc(global_time,0.d0,ps1,temp%startwbc,temp%nwbc)
822  end if
823  !!first step end
824 
825  evenstep=.true.
826 
827  tmpps2=>ps1
828 
829  do j=2,temp%s
830  bj(j)=dble(j**2+j-2)/dble(2*j*(j+1))
831  cmu=dble(2*j-1)/dble(j)*bj(j)/bj(j-1)
832  cmut=omega1*cmu
833  cnu=dble(1-j)/dble(j)*bj(j)/bj(j-2)
834  cnut=(bj(j-1)-1.d0)*cmut
835  one_mu_nu=1.d0-cmu-cnu
836  if(evenstep) then
837  tmpps1=>ps1
838  tmpps2=>ps2
839  else
840  tmpps1=>ps2
841  tmpps2=>ps1
842  end if
843 
844  dtj = cmut*my_dt
845  if(stagger_grid) then
846  !$OMP PARALLEL DO PRIVATE(igrid)
847  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
848  ! maybe the following global variables are needed in set_sources
849  ! next few lines ensure correct usage of routines like divvector etc
850  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
851  block=>ps(igrid)
852  ! end maybe the following global variables are needed in set_sources
853  call temp%sts_set_sources(ixg^ll,ixm^ll,tmpps1(igrid)%w,ps(igrid)%x,ps4(igrid)%w,fix_conserve_at_step,dtj,igrid,temp%nflux)
854  if(temp%nflux>ndir) then
855  tmpps2(igrid)%w(ixm^t,temp%startVar)=cmu*tmpps1(igrid)%w(ixm^t,temp%startVar)+&
856  cnu*tmpps2(igrid)%w(ixm^t,temp%startVar)+one_mu_nu*ps(igrid)%w(ixm^t,temp%startVar)+&
857  dtj*ps4(igrid)%w(ixm^t,temp%startVar)+cnut*ps3(igrid)%w(ixm^t,temp%startVar)
858  end if
859  tmpps2(igrid)%ws(ixc^s,1:nws)=cmu*tmpps1(igrid)%ws(ixc^s,1:nws)+&
860  cnu*tmpps2(igrid)%ws(ixc^s,1:nws)+one_mu_nu*ps(igrid)%ws(ixc^s,1:nws)+&
861  dtj*ps4(igrid)%w(ixc^s,iw_mag(1:nws))+cnut*ps3(igrid)%w(ixc^s,iw_mag(1:nws))
862  call phys_face_to_center(ixm^ll,tmpps2(igrid))
863  end do
864  !$OMP END PARALLEL DO
865  else
866  !$OMP PARALLEL DO PRIVATE(igrid)
867  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
868  ! maybe the following global variables are needed in set_sources
869  ! next few lines ensure correct usage of routines like divvector etc
870  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
871  block=>ps(igrid)
872  ! end maybe the following global variables are needed in set_sources
873  call temp%sts_set_sources(ixg^ll,ixm^ll,tmpps1(igrid)%w,ps(igrid)%x,ps4(igrid)%w,fix_conserve_at_step,dtj,igrid,temp%nflux)
874  tmpps2(igrid)%w(ixm^t,temp%startVar:temp%endVar)=cmu*tmpps1(igrid)%w(ixm^t,temp%startVar:temp%endVar)+&
875  cnu*tmpps2(igrid)%w(ixm^t,temp%startVar:temp%endVar)+one_mu_nu*ps(igrid)%w(ixm^t,temp%startVar:temp%endVar)+&
876  dtj*ps4(igrid)%w(ixm^t,temp%startVar:temp%endVar)+cnut*ps3(igrid)%w(ixm^t,temp%startVar:temp%endVar)
877  end do
878  !$OMP END PARALLEL DO
879  end if
880  if(fix_conserve_at_step) then
881  call recvflux(1,ndim)
882  call sendflux(1,ndim)
883  call fix_conserve(tmpps2,1,ndim,temp%startVar,temp%nflux)
884  if(stagger_grid) then
885  call fix_edges(tmpps2,1,ndim)
886  ! fill the cell-center values from the updated staggered variables
887  !$OMP PARALLEL DO PRIVATE(igrid)
888  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
889  call phys_face_to_center(ixg^ll,tmpps2(igrid))
890  end do
891  !$OMP END PARALLEL DO
892  end if
893  end if
894  if(associated(temp%sts_handle_errors)) then
895  !$OMP PARALLEL DO PRIVATE(igrid)
896  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
897  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
898  block=>ps(igrid)
899  call temp%sts_handle_errors(tmpps2(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,j)
900  end do
901  !$OMP END PARALLEL DO
902  end if
903  if(temp%nflux>temp%nwbc.and.temp%s==j) then
904  ! include the changed no-need-ghost-update variables in the last getbc
905  type_send_srl=>temp%type_send_srl_sts_2
906  type_recv_srl=>temp%type_recv_srl_sts_2
907  type_send_r=>temp%type_send_r_sts_2
908  type_recv_r=>temp%type_recv_r_sts_2
909  type_send_p=>temp%type_send_p_sts_2
910  type_recv_p=>temp%type_recv_p_sts_2
911  call getbc(global_time,0.d0,tmpps2,temp%startVar,temp%nflux)
912  else
913  call getbc(global_time,0.d0,tmpps2,temp%startwbc,temp%nwbc)
914  end if
915  evenstep=.not.evenstep
916  end do
917 
918  if(associated(temp%sts_after_last_cycle)) then
919  !$OMP PARALLEL DO PRIVATE(igrid)
920  do iigrid=1,igridstail; igrid=igrids(iigrid);
921  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
922  block=>ps(igrid)
923  ps(igrid)%w(ixg^t,temp%startVar:temp%endVar)=tmpps2(igrid)%w(ixg^t,temp%startVar:temp%endVar)
924  call temp%sts_after_last_cycle(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
925  end do
926  !$OMP END PARALLEL DO
927  phys_total_energy=total_energy_flag
928  prolongprimitive = prolong_flag
929  coarsenprimitive = coarsen_flag
930  else
931  if(stagger_grid) then
932  !$OMP PARALLEL DO PRIVATE(igrid)
933  do iigrid=1,igridstail; igrid=igrids(iigrid);
934  ps(igrid)%w(ixg^t,temp%startVar:temp%endVar)=tmpps2(igrid)%w(ixg^t,temp%startVar:temp%endVar)
935  ps(igrid)%ws=tmpps2(igrid)%ws
936  end do
937  !$OMP END PARALLEL DO
938  else
939  !$OMP PARALLEL DO PRIVATE(igrid)
940  do iigrid=1,igridstail; igrid=igrids(iigrid);
941  ps(igrid)%w(ixg^t,temp%startVar:temp%endVar)=tmpps2(igrid)%w(ixg^t,temp%startVar:temp%endVar)
942  end do
943  !$OMP END PARALLEL DO
944  end if
945  end if
946 
947  deallocate(bj)
948 
949  if(.not.temp%evolve_magnetic_field) then
950  ! restore stagger_grid value
951  stagger_grid=stagger_flag
952  end if
953 
954  temp=>temp%next
955  end do
956 
957  if(associated(head_sts_terms)) then
958  ! point bc mpi data type back to full type for (M)HD
965  end if
966 
967  bcphys=.true.
968 
969  if(phys_partial_ionization) then
970  ! update temperature variable in w
971  !$OMP PARALLEL DO PRIVATE(igrid)
972  do iigrid=1,igridstail; igrid=igrids(iigrid);
973  call phys_update_temperature(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%w,ps(igrid)%x)
974  end do
975  !$OMP END PARALLEL DO
976  end if
977 
978  end subroutine sts_add_source1
979 
980 end module mod_supertimestepping
interface for the function which gets the timestep in dtnew in the derived type
for the subroutines in this module, which do not depend on the term, but on the parameter sts_method ...
interface for the subroutines before_first_cycle and after_last_cycle in the derived type
interface for error handling subroutine in the derived type
subroutine, public alloc_state(igrid, s, ixGL, ixGextL, alloc_once_for_ps)
allocate memory to physical state of igrid node
Module for physical and numeric constants.
Definition: mod_constants.t:2
double precision, parameter dpi
Pi.
Definition: mod_constants.t:25
Module for flux conservation near refinement boundaries.
subroutine, public recvflux(idimLIM)
subroutine, public fix_edges(psuse, idimLIM)
subroutine, public sendflux(idimLIM)
subroutine, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
update ghost cells of all blocks including physical boundaries
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_f
integer, dimension(^nd, 0:3) l
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_f
integer, dimension(:^d &,:^d &), pointer type_recv_srl
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(:^d &,:^d &), pointer type_send_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_f
integer, dimension(:^d &,:^d &), pointer type_recv_r
integer, dimension(:^d &,:^d &), pointer type_send_r
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_f
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_f
integer, dimension(:^d &,:^d &), pointer type_recv_p
integer, dimension(:^d &,:^d &), pointer type_send_p
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
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.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision, dimension(:,:), allocatable dx
double precision, dimension(ndim) dxlevel
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
logical phys_partial_ionization
Solve partially ionized one-fluid plasma.
Definition: mod_physics.t:36
logical phys_total_energy
Solve total energy equation or not.
Definition: mod_physics.t:30
procedure(sub_update_temperature), pointer phys_update_temperature
Definition: mod_physics.t:92
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:83
Module for handling problematic values in simulations, such as negative pressures.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
integer, public sourcetype_sts
subroutine sts_add_source1(my_dt)
Iterates all the terms implemented with STS and adds the sources STS method 1 implementation.
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
integer function sts_get_ncycles1(dt, dtnew, dt_modified)
method used to set the number of cycles for the STS1 method
pure logical function, public is_sts_initialized()
integer, parameter, public sourcetype_sts_prior
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 ...
double precision function sum_chev(nu, N, limMax)
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
integer, parameter, public sourcetype_sts_split
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC....
pure double precision function chev(j, nu, N)
subroutine, public sts_init()
Initialize sts module.
integer, parameter, public sourcetype_sts_after
procedure(subr3), pointer, public sts_add_source