MPI-AMRVAC  3.0
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.9d0
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.5
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. 1d0) then
310  is=1
311  else
312  is=ceiling((dsqrt(9.d0+16.d0*ss)-1.d0)/2.d0)
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  end subroutine sts_add_source2
623 
624  !> Iterates all the terms implemented with STS and adds the sources
625  !> STS method 1 implementation
626  subroutine sts_add_source1(my_dt)
627  ! Meyer 2012 MNRAS 422,2102
630  use mod_fix_conserve
631  use mod_physics
632 
633  double precision, intent(in) :: my_dt
634  double precision :: dtj
635  double precision :: omega1,cmu,cmut,cnu,cnut,one_mu_nu
636  double precision, allocatable :: bj(:)
637  integer:: iigrid, igrid, j, ixC^L, ixGext^L
638  logical :: evenstep, stagger_flag=.false., prolong_flag=.false., coarsen_flag=.false., total_energy_flag=.true.
639  type(sts_term), pointer :: temp
640  type(state), dimension(:), pointer :: tmpPs1, tmpPs2
641 
642  ! do not fill physical boundary conditions
643  bcphys=.false.
644 
645  fix_conserve_at_step = time_advance .and. levmax>levmin
646 
647  temp => head_sts_terms
648  do while(associated(temp))
649 
650  if(.not.temp%evolve_magnetic_field) then
651  ! not do fix conserve and getbc for staggered values
652  stagger_flag=stagger_grid
653  stagger_grid=.false.
654  else if(stagger_grid) then
655  ixcmax^d=ixmhi^d;
656  ixcmin^d=ixmlo^d-1;
657  end if
658 
659  call init_comm_fix_conserve(1,ndim,temp%nflux)
660 
661  if(associated(temp%sts_before_first_cycle)) then
662  prolong_flag = prolongprimitive
663  coarsen_flag = coarsenprimitive
664  prolongprimitive = .false.
665  coarsenprimitive = .false.
666  total_energy_flag=phys_total_energy
667  phys_total_energy=.false.
668  !$OMP PARALLEL DO PRIVATE(igrid)
669  do iigrid=1,igridstail; igrid=igrids(iigrid);
670  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
671  block=>ps(igrid)
672  call temp%sts_before_first_cycle(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
673  if(.not. allocated(ps2(igrid)%w)) allocate(ps2(igrid)%w(ixg^t,1:nw))
674  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
675  if(.not. allocated(ps4(igrid)%w)) allocate(ps4(igrid)%w(ixg^t,1:nw))
676  ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
677  ps2(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
678  end do
679  !$OMP END PARALLEL DO
680  else
681  if(stagger_grid) then
682  ixgext^l=ixg^ll^ladd1;
683  !$OMP PARALLEL DO PRIVATE(igrid)
684  do iigrid=1,igridstail; igrid=igrids(iigrid);
685  if(.not. allocated(ps2(igrid)%w)) then
686  call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
687  end if
688  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
689  if(.not. allocated(ps4(igrid)%w)) allocate(ps4(igrid)%w(ixg^t,1:nw))
690  ps1(igrid)%w=ps(igrid)%w
691  ps2(igrid)%w=ps(igrid)%w
692  ps1(igrid)%ws=ps(igrid)%ws
693  ps2(igrid)%ws=ps(igrid)%ws
694  end do
695  !$OMP END PARALLEL DO
696  else
697  !$OMP PARALLEL DO PRIVATE(igrid)
698  do iigrid=1,igridstail; igrid=igrids(iigrid);
699  if(.not. allocated(ps2(igrid)%w)) allocate(ps2(igrid)%w(ixg^t,1:nw))
700  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
701  if(.not. allocated(ps4(igrid)%w)) allocate(ps4(igrid)%w(ixg^t,1:nw))
702  ps1(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
703  ps2(igrid)%w(ixg^t,1:nw)=ps(igrid)%w(ixg^t,1:nw)
704  end do
705  !$OMP END PARALLEL DO
706  end if
707  end if
708 
709  allocate(bj(0:temp%s))
710  bj(0)=1.d0/3.d0
711  bj(1)=bj(0)
712  if(temp%s>1) then
713  omega1=4.d0/dble(temp%s**2+temp%s-2)
714  cmut=omega1/3.d0
715  else
716  omega1=0.d0
717  cmut=1.d0
718  end if
719 
720  type_send_srl=>temp%type_send_srl_sts_1
721  type_recv_srl=>temp%type_recv_srl_sts_1
722  type_send_r=>temp%type_send_r_sts_1
723  type_recv_r=>temp%type_recv_r_sts_1
724  type_send_p=>temp%type_send_p_sts_1
725  type_recv_p=>temp%type_recv_p_sts_1
726 
727  if(.not. temp%types_initialized) then
728  call create_bc_mpi_datatype(temp%startwbc,temp%nwbc)
729  if(temp%nflux>temp%nwbc) then
730  ! prepare types for the changed no-need-ghost-update variables in the last getbc
731  type_send_srl=>temp%type_send_srl_sts_2
732  type_recv_srl=>temp%type_recv_srl_sts_2
733  type_send_r=>temp%type_send_r_sts_2
734  type_recv_r=>temp%type_recv_r_sts_2
735  type_send_p=>temp%type_send_p_sts_2
736  type_recv_p=>temp%type_recv_p_sts_2
737  call create_bc_mpi_datatype(temp%startVar,temp%nflux)
738  type_send_srl=>temp%type_send_srl_sts_1
739  type_recv_srl=>temp%type_recv_srl_sts_1
740  type_send_r=>temp%type_send_r_sts_1
741  type_recv_r=>temp%type_recv_r_sts_1
742  type_send_p=>temp%type_send_p_sts_1
743  type_recv_p=>temp%type_recv_p_sts_1
744  end if
745  temp%types_initialized = .true.
746  end if
747  dtj = cmut*my_dt
748  if(stagger_grid) then
749  !$OMP PARALLEL DO PRIVATE(igrid)
750  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
751  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
752  block=>ps(igrid)
753  ps4(igrid)%w=zero
754  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)
755  !!!eq solved: dU/dt = S, ps3 is stored S^n
756  ps3(igrid)%w(ixc^s,temp%startVar:temp%endVar) = my_dt * ps4(igrid)%w(ixc^s,temp%startVar:temp%endVar)
757  if(temp%nflux>ndir) then
758  ps1(igrid)%w(ixm^t,temp%startVar) = ps1(igrid)%w(ixm^t,temp%startVar) + cmut * ps3(igrid)%w(ixm^t,temp%startVar)
759  end if
760  ps1(igrid)%ws(ixc^s,1:nws) = ps1(igrid)%ws(ixc^s,1:nws) + cmut * ps3(igrid)%w(ixc^s,iw_mag(1:nws))
761  call phys_face_to_center(ixm^ll,ps1(igrid))
762  end do
763  !$OMP END PARALLEL DO
764  else
765  !$OMP PARALLEL DO PRIVATE(igrid)
766  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
767  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
768  block=>ps(igrid)
769  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)
770  !!!eq solved: dU/dt = S, ps3 is stored S^n
771  ps3(igrid)%w(ixm^t,temp%startVar:temp%endVar) = my_dt * ps4(igrid)%w(ixm^t,temp%startVar:temp%endVar)
772  ps1(igrid)%w(ixm^t,temp%startVar:temp%endVar) = ps1(igrid)%w(ixm^t,temp%startVar:temp%endVar) + &
773  cmut * ps3(igrid)%w(ixm^t,temp%startVar:temp%endVar)
774  end do
775  !$OMP END PARALLEL DO
776  end if
777  if(fix_conserve_at_step) then
778  call recvflux(1,ndim)
779  call sendflux(1,ndim)
780  call fix_conserve(ps1,1,ndim,temp%startVar,temp%nflux)
781  if(stagger_grid) then
782  call fix_edges(ps1,1,ndim)
783  ! fill the cell-center values from the updated staggered variables
784  !$OMP PARALLEL DO PRIVATE(igrid)
785  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
786  call phys_face_to_center(ixg^ll,ps1(igrid))
787  end do
788  !$OMP END PARALLEL DO
789  end if
790  end if
791  ! fix conservation of AMR grid by replacing flux from finer neighbors
792  if(associated(temp%sts_handle_errors)) then
793  !$OMP PARALLEL DO PRIVATE(igrid)
794  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
795  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
796  block=>ps(igrid)
797  call temp%sts_handle_errors(ps1(igrid)%w,ps1(igrid)%x,ixg^ll,ixm^ll,1)
798  end do
799  !$OMP END PARALLEL DO
800  end if
801  if(temp%nflux>temp%nwbc.and.temp%s==1) then
802  ! include the changed no-need-ghost-update variables in the last getbc
803  type_send_srl=>temp%type_send_srl_sts_2
804  type_recv_srl=>temp%type_recv_srl_sts_2
805  type_send_r=>temp%type_send_r_sts_2
806  type_recv_r=>temp%type_recv_r_sts_2
807  type_send_p=>temp%type_send_p_sts_2
808  type_recv_p=>temp%type_recv_p_sts_2
809  call getbc(global_time,0.d0,ps1,temp%startVar,temp%nflux)
810  else
811  call getbc(global_time,0.d0,ps1,temp%startwbc,temp%nwbc)
812  end if
813  !!first step end
814 
815  evenstep=.true.
816 
817  tmpps2=>ps1
818 
819  do j=2,temp%s
820  bj(j)=dble(j**2+j-2)/dble(2*j*(j+1))
821  cmu=dble(2*j-1)/dble(j)*bj(j)/bj(j-1)
822  cmut=omega1*cmu
823  cnu=dble(1-j)/dble(j)*bj(j)/bj(j-2)
824  cnut=(bj(j-1)-1.d0)*cmut
825  one_mu_nu=1.d0-cmu-cnu
826  if(evenstep) then
827  tmpps1=>ps1
828  tmpps2=>ps2
829  else
830  tmpps1=>ps2
831  tmpps2=>ps1
832  end if
833 
834  dtj = cmut*my_dt
835  if(stagger_grid) then
836  !$OMP PARALLEL DO PRIVATE(igrid)
837  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
838  ! maybe the following global variables are needed in set_sources
839  ! next few lines ensure correct usage of routines like divvector etc
840  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
841  block=>ps(igrid)
842  ! end maybe the following global variables are needed in set_sources
843  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)
844  if(temp%nflux>ndir) then
845  tmpps2(igrid)%w(ixm^t,temp%startVar)=cmu*tmpps1(igrid)%w(ixm^t,temp%startVar)+&
846  cnu*tmpps2(igrid)%w(ixm^t,temp%startVar)+one_mu_nu*ps(igrid)%w(ixm^t,temp%startVar)+&
847  dtj*ps4(igrid)%w(ixm^t,temp%startVar)+cnut*ps3(igrid)%w(ixm^t,temp%startVar)
848  end if
849  tmpps2(igrid)%ws(ixc^s,1:nws)=cmu*tmpps1(igrid)%ws(ixc^s,1:nws)+&
850  cnu*tmpps2(igrid)%ws(ixc^s,1:nws)+one_mu_nu*ps(igrid)%ws(ixc^s,1:nws)+&
851  dtj*ps4(igrid)%w(ixc^s,iw_mag(1:nws))+cnut*ps3(igrid)%w(ixc^s,iw_mag(1:nws))
852  call phys_face_to_center(ixm^ll,tmpps2(igrid))
853  end do
854  !$OMP END PARALLEL DO
855  else
856  !$OMP PARALLEL DO PRIVATE(igrid)
857  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
858  ! maybe the following global variables are needed in set_sources
859  ! next few lines ensure correct usage of routines like divvector etc
860  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
861  block=>ps(igrid)
862  ! end maybe the following global variables are needed in set_sources
863  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)
864  tmpps2(igrid)%w(ixm^t,temp%startVar:temp%endVar)=cmu*tmpps1(igrid)%w(ixm^t,temp%startVar:temp%endVar)+&
865  cnu*tmpps2(igrid)%w(ixm^t,temp%startVar:temp%endVar)+one_mu_nu*ps(igrid)%w(ixm^t,temp%startVar:temp%endVar)+&
866  dtj*ps4(igrid)%w(ixm^t,temp%startVar:temp%endVar)+cnut*ps3(igrid)%w(ixm^t,temp%startVar:temp%endVar)
867  end do
868  !$OMP END PARALLEL DO
869  end if
870  if(fix_conserve_at_step) then
871  call recvflux(1,ndim)
872  call sendflux(1,ndim)
873  call fix_conserve(tmpps2,1,ndim,temp%startVar,temp%nflux)
874  if(stagger_grid) then
875  call fix_edges(tmpps2,1,ndim)
876  ! fill the cell-center values from the updated staggered variables
877  !$OMP PARALLEL DO PRIVATE(igrid)
878  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
879  call phys_face_to_center(ixg^ll,tmpps2(igrid))
880  end do
881  !$OMP END PARALLEL DO
882  end if
883  end if
884  if(associated(temp%sts_handle_errors)) then
885  !$OMP PARALLEL DO PRIVATE(igrid)
886  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
887  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
888  block=>ps(igrid)
889  call temp%sts_handle_errors(tmpps2(igrid)%w,ps(igrid)%x,ixg^ll,ixm^ll,j)
890  end do
891  !$OMP END PARALLEL DO
892  end if
893  if(temp%nflux>temp%nwbc.and.temp%s==j) then
894  ! include the changed no-need-ghost-update variables in the last getbc
895  type_send_srl=>temp%type_send_srl_sts_2
896  type_recv_srl=>temp%type_recv_srl_sts_2
897  type_send_r=>temp%type_send_r_sts_2
898  type_recv_r=>temp%type_recv_r_sts_2
899  type_send_p=>temp%type_send_p_sts_2
900  type_recv_p=>temp%type_recv_p_sts_2
901  call getbc(global_time,0.d0,tmpps2,temp%startVar,temp%nflux)
902  else
903  call getbc(global_time,0.d0,tmpps2,temp%startwbc,temp%nwbc)
904  end if
905  evenstep=.not.evenstep
906  end do
907 
908  if(associated(temp%sts_after_last_cycle)) then
909  !$OMP PARALLEL DO PRIVATE(igrid)
910  do iigrid=1,igridstail; igrid=igrids(iigrid);
911  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
912  block=>ps(igrid)
913  ps(igrid)%w(ixg^t,temp%startVar:temp%endVar)=tmpps2(igrid)%w(ixg^t,temp%startVar:temp%endVar)
914  call temp%sts_after_last_cycle(ixg^ll,ixg^ll,ps(igrid)%w,ps(igrid)%x)
915  end do
916  !$OMP END PARALLEL DO
917  phys_total_energy=total_energy_flag
918  prolongprimitive = prolong_flag
919  coarsenprimitive = coarsen_flag
920  else
921  if(stagger_grid) then
922  !$OMP PARALLEL DO PRIVATE(igrid)
923  do iigrid=1,igridstail; igrid=igrids(iigrid);
924  ps(igrid)%w(ixg^t,temp%startVar:temp%endVar)=tmpps2(igrid)%w(ixg^t,temp%startVar:temp%endVar)
925  ps(igrid)%ws=tmpps2(igrid)%ws
926  end do
927  !$OMP END PARALLEL DO
928  else
929  !$OMP PARALLEL DO PRIVATE(igrid)
930  do iigrid=1,igridstail; igrid=igrids(iigrid);
931  ps(igrid)%w(ixg^t,temp%startVar:temp%endVar)=tmpps2(igrid)%w(ixg^t,temp%startVar:temp%endVar)
932  end do
933  !$OMP END PARALLEL DO
934  end if
935  endif
936 
937  deallocate(bj)
938 
939  if(.not.temp%evolve_magnetic_field) then
940  ! restore stagger_grid value
941  stagger_grid=stagger_flag
942  end if
943 
944  temp=>temp%next
945  end do
946 
947  if(associated(head_sts_terms)) then
948  ! point bc mpi data type back to full type for (M)HD
955  end if
956 
957  bcphys=.true.
958 
959  end subroutine sts_add_source1
960 
961 end module mod_supertimestepping
subroutine alloc_state(igrid, s, ixGL, ixGextL, alloc_once_for_ps)
allocate memory to physical state of igrid node
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
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
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 (within a block with 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_total_energy
Solve total energy equation or not.
Definition: mod_physics.t:33
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:84
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