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