MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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
35 !The following is used only for method 2, not input parameter TODO check if we want as input parameter
36 double precision,parameter :: nu_sts = 0.5d0
37 !> the maximum number of subcycles
38 integer :: sts_ncycles=1000
39 integer :: sts_method = 1
40 integer, parameter :: sourcetype_sts_prior =0
41 integer, parameter :: sourcetype_sts_after =1
42 integer, parameter :: sourcetype_sts_split =2
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 double precision :: dt_expl
107 integer, public :: s
108
109 !>types used for send/recv ghosts, see mod_ghostcells_update
110 integer, dimension(-1:1^D&) :: type_send_srl_sts_1, type_recv_srl_sts_1
111 integer, dimension(-1:1^D&) :: type_send_r_sts_1
112 integer, dimension( 0:3^D&) :: type_recv_r_sts_1
113 integer, dimension( 0:3^D&) :: type_recv_p_sts_1, type_send_p_sts_1
114
115 integer, dimension(-1:1^D&) :: type_send_srl_sts_2, type_recv_srl_sts_2
116 integer, dimension(-1:1^D&) :: type_send_r_sts_2
117 integer, dimension( 0:3^D&) :: type_recv_r_sts_2
118 integer, dimension( 0:3^D&) :: type_recv_p_sts_2, type_send_p_sts_2
119
120 integer :: startvar
121 integer :: endvar
122 !> number of flux involved in STS update
123 integer :: nflux
124 integer :: startwbc
125 integer :: nwbc
126 logical :: types_initialized
127 logical :: evolve_magnetic_field
128 procedure(subr1), pointer, nopass :: sts_set_sources
129 procedure(subr2), pointer, nopass :: sts_getdt
130 procedure(subr5), pointer, nopass :: sts_before_first_cycle, sts_after_last_cycle
131 procedure(subr_e), pointer, nopass :: sts_handle_errors
132 type(sts_term), pointer :: next
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
142contains
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
154 sts_add_source => sts_add_source1
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)
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)
186111 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)
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)
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 double precision :: dtnew,dtmin_mype
353 double precision :: dx^d, ss
354 integer:: iigrid, igrid, ncycles
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
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
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
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
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
980end module mod_supertimestepping
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.
double precision, parameter dpi
Pi.
Module for flux conservation near refinement boundaries.
subroutine, public init_comm_fix_conserve(idimlim, nwfluxin)
subroutine, public fix_edges(psuse, idimlim)
subroutine, public recvflux(idimlim)
subroutine, public sendflux(idimlim)
subroutine, public fix_conserve(psb, idimlim, nw0, 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( :^d &), pointer type_recv_r
integer, dimension(^nd, 0:3) l
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_send_p_f
integer, dimension( :^d &), pointer type_send_p
integer, dimension( :^d &), pointer type_send_srl
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(0:3^d &), target type_recv_r_f
integer, dimension(-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &), target type_send_r_f
integer, dimension(-1:1^d &), target type_send_srl_f
integer, dimension( :^d &), pointer type_send_r
integer, dimension( :^d &), pointer type_recv_p
integer, dimension(0:3^d &), target type_recv_p_f
integer, dimension( :^d &), pointer type_recv_srl
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.
logical coarsenprimitive
coarsen primitive variables in level-jump ghost cells
double precision, 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.
logical time_advance
do time evolving
logical prolongprimitive
prolongate primitive variables in level-jump ghost cells
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision, dimension(:,:), allocatable dx
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
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:44
logical phys_total_energy
Solve total energy equation or not.
Definition mod_physics.t:38
procedure(sub_update_temperature), pointer phys_update_temperature
Definition mod_physics.t:90
procedure(sub_face_to_center), pointer phys_face_to_center
Definition mod_physics.t:81
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...
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 ...
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...
type(sts_term), pointer head_sts_terms
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....
subroutine, public sts_init()
Initialize sts module.
integer, parameter, public sourcetype_sts_after
procedure(subr3), pointer, public sts_add_source