MPI-AMRVAC  3.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_ghostcells_update.t
Go to the documentation of this file.
1 !> update ghost cells of all blocks including physical boundaries
3 
4  implicit none
5  ! Special buffer for pole copy
6  type wbuffer
7  double precision, dimension(:^D&,:), allocatable :: w
8  end type wbuffer
9 
10  ! A switch of update physical boundary or not
11  logical, public :: bcphys=.true.
12 
13  integer :: ixm^l, ixcog^l, ixcom^l, ixcogs^l
14 
15  ! The number of interleaving sending buffers for ghost cells
16  integer, parameter :: npwbuf=2
17 
18  ! The first index goes from -1:2, where -1 is used when a block touches the
19  ! lower boundary, 1 when a block touches an upper boundary, and 0 a situation
20  ! away from boundary conditions, 2 when a block touched both lower and upper
21  ! boundary
22 
23  ! index ranges to send (S) to sibling blocks, receive (R) from sibling blocks
24  integer, dimension(-1:2,-1:1) :: ixs_srl_^l, ixr_srl_^l
25 
26  ! index ranges of staggered variables to send (S) to sibling blocks, receive (R) from sibling blocks
27  integer, dimension(^ND,-1:1) :: ixs_srl_stg_^l, ixr_srl_stg_^l
28 
29  ! index ranges to send (S) restricted (r) ghost cells to coarser blocks
30  integer, dimension(-1:1,-1:1) :: ixs_r_^l
31 
32  ! index ranges of staggered variables to send (S) restricted (r) ghost cells to coarser blocks
33  integer, dimension(^ND,-1:1) :: ixs_r_stg_^l
34 
35  ! index ranges to receive restriced ghost cells from finer blocks
36  integer, dimension(-1:1, 0:3) :: ixr_r_^l
37 
38  ! index ranges of staggered variables to receive restriced ghost cells from finer blocks
39  integer, dimension(^ND,0:3) :: ixr_r_stg_^l
40 
41  ! send prolongated (p) ghost cells to finer blocks, receive prolongated from coarser blocks
42  integer, dimension(-1:1, 0:3) :: ixs_p_^l, ixr_p_^l
43 
44  ! send prolongated (p) staggered ghost cells to finer blocks, receive prolongated from coarser blocks
45  integer, dimension(^ND,0:3) :: ixs_p_stg_^l, ixr_p_stg_^l
46 
47  ! number of MPI receive-send pairs, srl: same refinement level; r: restrict; p: prolong
49 
50  ! record index position of buffer arrays
52 
53  ! count of times of send and receive
55 
56  ! count of times of send and receive for cell center ghost cells
57  integer :: isend_c, irecv_c
58 
59  ! tag of MPI send and recv
60  integer, private :: itag
61 
62  ! total sizes = cell-center normal flux + stagger-grid flux of send and receive
63  integer, dimension(-1:1^D&) :: sizes_srl_send_total, sizes_srl_recv_total
64 
65  ! sizes of buffer arrays for center-grid variable for siblings and restrict
66  integer, dimension(:), allocatable :: recvrequest_c_sr, sendrequest_c_sr
67  integer, dimension(:,:), allocatable :: recvstatus_c_sr, sendstatus_c_sr
68 
69  ! sizes of buffer arrays for center-grid variable for prolongation
70  integer, dimension(:), allocatable :: recvrequest_c_p, sendrequest_c_p
71  integer, dimension(:,:), allocatable :: recvstatus_c_p, sendstatus_c_p
72 
73  ! sizes of buffer arrays for stagger-grid variable
74  integer, dimension(^ND,-1:1^D&) :: sizes_srl_send_stg, sizes_srl_recv_stg
75 
76  integer, dimension(:), allocatable :: recvrequest_srl, sendrequest_srl
77  integer, dimension(:,:), allocatable :: recvstatus_srl, sendstatus_srl
78 
79  ! buffer arrays for send and receive of siblings, allocate in build_connectivity
80  double precision, dimension(:), allocatable :: recvbuffer_srl, sendbuffer_srl
81 
82  integer, dimension(:), allocatable :: recvrequest_r, sendrequest_r
83  integer, dimension(:,:), allocatable :: recvstatus_r, sendstatus_r
84 
85  ! buffer arrays for send and receive in restriction
86  double precision, dimension(:), allocatable :: recvbuffer_r, sendbuffer_r
87 
88  integer, dimension(:), allocatable :: recvrequest_p, sendrequest_p
89  integer, dimension(:,:), allocatable :: recvstatus_p, sendstatus_p
90 
91  ! buffer arrays for send and receive in prolongation
92  double precision, dimension(:), allocatable :: recvbuffer_p, sendbuffer_p
93 
94  ! sizes to allocate buffer arrays for send and receive for restriction
95  integer, dimension(-1:1^D&) :: sizes_r_send_total
96  integer, dimension(0:3^D&) :: sizes_r_recv_total
97  integer, dimension(^ND,-1:1^D&) :: sizes_r_send_stg
98  integer, dimension(^ND,0:3^D&) :: sizes_r_recv_stg
99 
100  ! sizes to allocate buffer arrays for send and receive for restriction
101  integer, dimension(0:3^D&) :: sizes_p_send_total, sizes_p_recv_total
102  integer, dimension(^ND,0:3^D&) :: sizes_p_send_stg, sizes_p_recv_stg
103 
104  ! There are two variants, _f indicates that all flux variables are filled,
105  ! whereas _p means that part of the variables is filled
106  ! Furthermore _r_ stands for restrict, _p_ for prolongation.
107  integer, dimension(-1:2^D&,-1:1^D&), target :: type_send_srl_f, type_recv_srl_f
108  integer, dimension(-1:1^D&,-1:1^D&), target :: type_send_r_f
109  integer, dimension(-1:1^D&, 0:3^D&), target :: type_recv_r_f, type_send_p_f, type_recv_p_f
110  integer, dimension(-1:2^D&,-1:1^D&), target :: type_send_srl_p1, type_recv_srl_p1
111  integer, dimension(-1:1^D&,-1:1^D&), target :: type_send_r_p1
112  integer, dimension(-1:1^D&, 0:3^D&), target :: type_recv_r_p1, type_send_p_p1, type_recv_p_p1
113  integer, dimension(-1:2^D&,-1:1^D&), target :: type_send_srl_p2, type_recv_srl_p2
114  integer, dimension(-1:1^D&,-1:1^D&), target :: type_send_r_p2
115  integer, dimension(-1:1^D&, 0:3^D&), target :: type_recv_r_p2, type_send_p_p2, type_recv_p_p2
116  integer, dimension(:^D&,:^D&), pointer :: type_send_srl, type_recv_srl, type_send_r
117  integer, dimension(:^D&,:^D&), pointer :: type_recv_r, type_send_p, type_recv_p
118 
119 contains
120 
121  subroutine init_bc()
124 
125  integer :: nghostcellsCo, interpolation_order
126  integer :: nx^D, nxCo^D, ixG^L, i^D, ic^D, inc^D, idir
127 
128  ixg^l=ixg^ll;
129  ixm^l=ixg^l^lsubnghostcells;
130  ixcogmin^d=1;
131  ixcogmax^d=(ixghi^d-2*nghostcells)/2+2*nghostcells;
132  ixcogsmin^d=0;
133  ixcogsmax^d=ixcogmax^d;
134 
135  ixcom^l=ixcog^l^lsubnghostcells;
136 
137  nx^d=ixmmax^d-ixmmin^d+1;
138  nxco^d=nx^d/2;
139 
140  if(ghost_copy) then
141  interpolation_order=1
142  else
143  interpolation_order=2
144  end if
145  nghostcellsco=int((nghostcells+1)/2)
146 
147  if (nghostcellsco+interpolation_order-1>nghostcells) then
148  call mpistop("interpolation order for prolongation in getbc too high")
149  end if
150 
151  ! (iib,i) index has following meanings: iib = 0 means it is not at any physical boundary
152  ! iib=-1 means it is at the minimum side of a physical boundary
153  ! iib= 1 means it is at the maximum side of a physical boundary
154  ! i=-1 means subregion prepared for the neighbor at its minimum side
155  ! i= 1 means subregion prepared for the neighbor at its maximum side
156  {
157  ixs_srl_min^d(:,-1)=ixmmin^d
158  ixs_srl_min^d(:, 0)=ixmmin^d
159  ixs_srl_min^d(:, 1)=ixmmax^d+1-nghostcells
160  ixs_srl_max^d(:,-1)=ixmmin^d-1+nghostcells
161  ixs_srl_max^d(:, 0)=ixmmax^d
162  ixs_srl_max^d(:, 1)=ixmmax^d
163 
164  ixr_srl_min^d(:,-1)=1
165  ixr_srl_min^d(:, 0)=ixmmin^d
166  ixr_srl_min^d(:, 1)=ixmmax^d+1
167  ixr_srl_max^d(:,-1)=nghostcells
168  ixr_srl_max^d(:, 0)=ixmmax^d
169  ixr_srl_max^d(:, 1)=ixgmax^d
170 
171  ixs_r_min^d(:,-1)=ixcommin^d
172  ixs_r_min^d(:, 0)=ixcommin^d
173  ixs_r_min^d(:, 1)=ixcommax^d+1-nghostcells
174  ixs_r_max^d(:,-1)=ixcommin^d-1+nghostcells
175  ixs_r_max^d(:, 0)=ixcommax^d
176  ixs_r_max^d(:, 1)=ixcommax^d
177 
178  ixr_r_min^d(:, 0)=1
179  ixr_r_min^d(:, 1)=ixmmin^d
180  ixr_r_min^d(:, 2)=ixmmin^d+nxco^d
181  ixr_r_min^d(:, 3)=ixmmax^d+1
182  ixr_r_max^d(:, 0)=nghostcells
183  ixr_r_max^d(:, 1)=ixmmin^d-1+nxco^d
184  ixr_r_max^d(:, 2)=ixmmax^d
185  ixr_r_max^d(:, 3)=ixgmax^d
186 
187  ixs_p_min^d(:, 0)=ixmmin^d-(interpolation_order-1)
188  ixs_p_min^d(:, 1)=ixmmin^d-(interpolation_order-1)
189  ixs_p_min^d(:, 2)=ixmmin^d+nxco^d-nghostcellsco-(interpolation_order-1)
190  ixs_p_min^d(:, 3)=ixmmax^d+1-nghostcellsco-(interpolation_order-1)
191  ixs_p_max^d(:, 0)=ixmmin^d-1+nghostcellsco+(interpolation_order-1)
192  ixs_p_max^d(:, 1)=ixmmin^d-1+nxco^d+nghostcellsco+(interpolation_order-1)
193  ixs_p_max^d(:, 2)=ixmmax^d+(interpolation_order-1)
194  ixs_p_max^d(:, 3)=ixmmax^d+(interpolation_order-1)
195 
196  if(.not.phys_req_diagonal) then
197  ! exclude ghost-cell region when diagonal cells are unknown
198  ixs_p_min^d(:, 0)=ixmmin^d
199  ixs_p_max^d(:, 3)=ixmmax^d
200  ixs_p_max^d(:, 1)=ixmmin^d-1+nxco^d+(interpolation_order-1)
201  ixs_p_min^d(:, 2)=ixmmin^d+nxco^d-(interpolation_order-1)
202  end if
203 
204  ixr_p_min^d(:, 0)=ixcommin^d-nghostcellsco-(interpolation_order-1)
205  ixr_p_min^d(:, 1)=ixcommin^d-(interpolation_order-1)
206  ixr_p_min^d(:, 2)=ixcommin^d-nghostcellsco-(interpolation_order-1)
207  ixr_p_min^d(:, 3)=ixcommax^d+1-(interpolation_order-1)
208  ixr_p_max^d(:, 0)=nghostcells+(interpolation_order-1)
209  ixr_p_max^d(:, 1)=ixcommax^d+nghostcellsco+(interpolation_order-1)
210  ixr_p_max^d(:, 2)=ixcommax^d+(interpolation_order-1)
211  ixr_p_max^d(:, 3)=ixcommax^d+nghostcellsco+(interpolation_order-1)
212 
213  if(.not.phys_req_diagonal) then
214  ixr_p_max^d(:, 0)=nghostcells
215  ixr_p_min^d(:, 3)=ixcommax^d+1
216  ixr_p_max^d(:, 1)=ixcommax^d+(interpolation_order-1)
217  ixr_p_min^d(:, 2)=ixcommin^d-(interpolation_order-1)
218  end if
219 
220  \}
221 
222  if (stagger_grid) then
223  allocate(pole_buf%ws(ixgs^t,nws))
224  ! Staggered (face-allocated) variables
225  do idir=1,ndim
226  { ixs_srl_stg_min^d(idir,-1)=ixmmin^d-kr(idir,^d)
227  ixs_srl_stg_max^d(idir,-1)=ixmmin^d-1+nghostcells
228  ixs_srl_stg_min^d(idir,0) =ixmmin^d-kr(idir,^d)
229  ixs_srl_stg_max^d(idir,0) =ixmmax^d
230  ixs_srl_stg_min^d(idir,1) =ixmmax^d-nghostcells+1-kr(idir,^d)
231  ixs_srl_stg_max^d(idir,1) =ixmmax^d
232 
233  ixr_srl_stg_min^d(idir,-1)=1-kr(idir,^d)
234  ixr_srl_stg_max^d(idir,-1)=nghostcells
235  ixr_srl_stg_min^d(idir,0) =ixmmin^d-kr(idir,^d)
236  ixr_srl_stg_max^d(idir,0) =ixmmax^d
237  ixr_srl_stg_min^d(idir,1) =ixmmax^d+1-kr(idir,^d)
238  ixr_srl_stg_max^d(idir,1) =ixgmax^d
239 
240  ixs_r_stg_min^d(idir,-1)=ixcommin^d-kr(idir,^d)
241  ixs_r_stg_max^d(idir,-1)=ixcommin^d-1+nghostcells
242  ixs_r_stg_min^d(idir,0) =ixcommin^d-kr(idir,^d)
243  ixs_r_stg_max^d(idir,0) =ixcommax^d
244  ixs_r_stg_min^d(idir,1) =ixcommax^d+1-nghostcells-kr(idir,^d)
245  ixs_r_stg_max^d(idir,1) =ixcommax^d
246 
247  ixr_r_stg_min^d(idir,0)=1-kr(idir,^d)
248  ixr_r_stg_max^d(idir,0)=nghostcells
249  ixr_r_stg_min^d(idir,1)=ixmmin^d-kr(idir,^d)
250  ixr_r_stg_max^d(idir,1)=ixmmin^d-1+nxco^d
251  ixr_r_stg_min^d(idir,2)=ixmmin^d+nxco^d-kr(idir,^d)
252  ixr_r_stg_max^d(idir,2)=ixmmax^d
253  ixr_r_stg_min^d(idir,3)=ixmmax^d+1-kr(idir,^d)
254  ixr_r_stg_max^d(idir,3)=ixgmax^d
255  \}
256  {if (idir==^d) then
257  ! Parallel components
258  {
259  ixs_p_stg_min^d(idir,0)=ixmmin^d-1 ! -1 to make redundant
260  ixs_p_stg_max^d(idir,0)=ixmmin^d-1+nghostcellsco
261  ixs_p_stg_min^d(idir,1)=ixmmin^d-1 ! -1 to make redundant
262  ixs_p_stg_max^d(idir,1)=ixmmin^d-1+nxco^d+nghostcellsco
263  ixs_p_stg_min^d(idir,2)=ixmmax^d-nxco^d-nghostcellsco
264  ixs_p_stg_max^d(idir,2)=ixmmax^d
265  ixs_p_stg_min^d(idir,3)=ixmmax^d-nghostcellsco
266  ixs_p_stg_max^d(idir,3)=ixmmax^d
267 
268  ixr_p_stg_min^d(idir,0)=ixcommin^d-1-nghostcellsco
269  ixr_p_stg_max^d(idir,0)=ixcommin^d-1
270  ixr_p_stg_min^d(idir,1)=ixcommin^d-1 ! -1 to make redundant
271  ixr_p_stg_max^d(idir,1)=ixcommax^d+nghostcellsco
272  ixr_p_stg_min^d(idir,2)=ixcommin^d-1-nghostcellsco
273  ixr_p_stg_max^d(idir,2)=ixcommax^d
274  ixr_p_stg_min^d(idir,3)=ixcommax^d+1-1 ! -1 to make redundant
275  ixr_p_stg_max^d(idir,3)=ixcommax^d+nghostcellsco
276  \}
277  else
278  {
279  ! Perpendicular component
280  ixs_p_stg_min^d(idir,0)=ixmmin^d
281  ixs_p_stg_max^d(idir,0)=ixmmin^d-1+nghostcellsco+(interpolation_order-1)
282  ixs_p_stg_min^d(idir,1)=ixmmin^d
283  ixs_p_stg_max^d(idir,1)=ixmmin^d-1+nxco^d+nghostcellsco+(interpolation_order-1)
284  ixs_p_stg_min^d(idir,2)=ixmmax^d+1-nxco^d-nghostcellsco-(interpolation_order-1)
285  ixs_p_stg_max^d(idir,2)=ixmmax^d
286  ixs_p_stg_min^d(idir,3)=ixmmax^d+1-nghostcellsco-(interpolation_order-1)
287  ixs_p_stg_max^d(idir,3)=ixmmax^d
288 
289  ixr_p_stg_min^d(idir,0)=ixcommin^d-nghostcellsco-(interpolation_order-1)
290  ixr_p_stg_max^d(idir,0)=ixcommin^d-1
291  ixr_p_stg_min^d(idir,1)=ixcommin^d
292  ixr_p_stg_max^d(idir,1)=ixcommax^d+nghostcellsco+(interpolation_order-1)
293  ixr_p_stg_min^d(idir,2)=ixcommin^d-nghostcellsco-(interpolation_order-1)
294  ixr_p_stg_max^d(idir,2)=ixcommax^d
295  ixr_p_stg_min^d(idir,3)=ixcommax^d+1
296  ixr_p_stg_max^d(idir,3)=ixcommax^d+nghostcellsco+(interpolation_order-1)
297  \}
298  end if
299  }
300  end do
301  ! calculate sizes for buffer arrays for siblings
302  {do i^db=-1,1\}
303  ! Staggered (face-allocated) variables
304  do idir=1,ndim
305  sizes_srl_send_stg(idir,i^d)={(ixs_srl_stg_max^d(idir,i^d)-ixs_srl_stg_min^d(idir,i^d)+1)|*}
306  sizes_srl_recv_stg(idir,i^d)={(ixr_srl_stg_max^d(idir,i^d)-ixr_srl_stg_min^d(idir,i^d)+1)|*}
307  sizes_r_send_stg(idir,i^d)={(ixs_r_stg_max^d(idir,i^d)-ixs_r_stg_min^d(idir,i^d)+1)|*}
308  end do
311  sizes_r_send_total(i^d)=sum(sizes_r_send_stg(:,i^d))
312  {end do\}
313 
314  {do i^db=0,3\}
315  ! Staggered (face-allocated) variables
316  do idir=1,ndim
317  sizes_r_recv_stg(idir,i^d)={(ixr_r_stg_max^d(idir,i^d)-ixr_r_stg_min^d(idir,i^d)+1)|*}
318  sizes_p_send_stg(idir,i^d)={(ixs_p_stg_max^d(idir,i^d)-ixs_p_stg_min^d(idir,i^d)+1)|*}
319  sizes_p_recv_stg(idir,i^d)={(ixr_p_stg_max^d(idir,i^d)-ixr_p_stg_min^d(idir,i^d)+1)|*}
320  end do
321  sizes_r_recv_total(i^d)=sum(sizes_r_recv_stg(:,i^d))
322  sizes_p_send_total(i^d)=sum(sizes_p_send_stg(:,i^d))
323  sizes_p_recv_total(i^d)=sum(sizes_p_recv_stg(:,i^d))
324  {end do\}
325  end if
326  if(.not.stagger_grid .or. physics_type=='mf') then
327  ! extend index range to physical boundary
328  {
329  ixs_srl_min^d(-1,0)=1
330  ixs_srl_min^d( 1,0)=ixmmin^d
331  ixs_srl_min^d( 2,0)=1
332  ixs_srl_max^d(-1,0)=ixmmax^d
333  ixs_srl_max^d( 1,0)=ixgmax^d
334  ixs_srl_max^d( 2,0)=ixgmax^d
335 
336  ixr_srl_min^d(-1,0)=1
337  ixr_srl_min^d( 1,0)=ixmmin^d
338  ixr_srl_min^d( 2,0)=1
339  ixr_srl_max^d(-1,0)=ixmmax^d
340  ixr_srl_max^d( 1,0)=ixgmax^d
341  ixr_srl_max^d( 2,0)=ixgmax^d
342 
343  ixs_r_min^d(-1,0)=1
344  ixs_r_min^d( 1,0)=ixcommin^d
345  ixs_r_max^d(-1,0)=ixcommax^d
346  ixs_r_max^d( 1,0)=ixcogmax^d
347 
348  ixr_r_min^d(-1,1)=1
349  ixr_r_max^d(-1,1)=ixmmin^d-1+nxco^d
350  ixr_r_min^d( 1,2)=ixmmin^d+nxco^d
351  ixr_r_max^d( 1,2)=ixgmax^d
352 
353  ixs_p_min^d(-1,1)=1
354  ixs_p_max^d( 1,2)=ixgmax^d
355 
356  ixr_p_min^d(-1,1)=1
357  ixr_p_max^d( 1,2)=ixcogmax^d
358  \}
359  end if
360 
361  end subroutine init_bc
362 
363  subroutine create_bc_mpi_datatype(nwstart,nwbc)
365 
366  integer, intent(in) :: nwstart, nwbc
367  integer :: i^D, ic^D, inc^D, iib^D
368 
369  {do i^db=-1,1\}
370  if (i^d==0|.and.) cycle
371  {do iib^db=-1,2\}
372  call get_bc_comm_type(type_send_srl(iib^d,i^d),ixs_srl_^l(iib^d,i^d),ixg^ll,nwstart,nwbc)
373  call get_bc_comm_type(type_recv_srl(iib^d,i^d),ixr_srl_^l(iib^d,i^d),ixg^ll,nwstart,nwbc)
374  if (iib^d==2|.or.) cycle
375  call get_bc_comm_type(type_send_r(iib^d,i^d), ixs_r_^l(iib^d,i^d),ixcog^l,nwstart,nwbc)
376  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
377  inc^db=2*i^db+ic^db\}
378  call get_bc_comm_type(type_recv_r(iib^d,inc^d),ixr_r_^l(iib^d,inc^d), ixg^ll,nwstart,nwbc)
379  call get_bc_comm_type(type_send_p(iib^d,inc^d),ixs_p_^l(iib^d,inc^d), ixg^ll,nwstart,nwbc)
380  call get_bc_comm_type(type_recv_p(iib^d,inc^d),ixr_p_^l(iib^d,inc^d),ixcog^l,nwstart,nwbc)
381  {end do\}
382  {end do\}
383  {end do\}
384 
385  end subroutine create_bc_mpi_datatype
386 
387  subroutine get_bc_comm_type(comm_type,ix^L,ixG^L,nwstart,nwbc)
389 
390  integer, intent(inout) :: comm_type
391  integer, intent(in) :: ix^L, ixG^L, nwstart, nwbc
392 
393  integer, dimension(ndim+1) :: fullsize, subsize, start
394 
395  ^d&fullsize(^d)=ixgmax^d;
396  fullsize(ndim+1)=nw
397  ^d&subsize(^d)=ixmax^d-ixmin^d+1;
398  subsize(ndim+1)=nwbc
399  ^d&start(^d)=ixmin^d-1;
400  start(ndim+1)=nwstart-1
401 
402  call mpi_type_create_subarray(ndim+1,fullsize,subsize,start,mpi_order_fortran, &
403  mpi_double_precision,comm_type,ierrmpi)
404  call mpi_type_commit(comm_type,ierrmpi)
405 
406  end subroutine get_bc_comm_type
407 
408  subroutine put_bc_comm_types()
410 
411  integer :: i^D, ic^D, inc^D, iib^D
412 
413  {do i^db=-1,1\}
414  if (i^d==0|.and.) cycle
415  {do iib^db=-1,2\}
416  call mpi_type_free(type_send_srl(iib^d,i^d),ierrmpi)
417  call mpi_type_free(type_recv_srl(iib^d,i^d),ierrmpi)
418  if (levmin==levmax) cycle
419  if (iib^d==2|.or.) cycle
420  call mpi_type_free(type_send_r(iib^d,i^d),ierrmpi)
421  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
422  inc^db=2*i^db+ic^db\}
423  call mpi_type_free(type_recv_r(iib^d,inc^d),ierrmpi)
424  call mpi_type_free(type_send_p(iib^d,inc^d),ierrmpi)
425  call mpi_type_free(type_recv_p(iib^d,inc^d),ierrmpi)
426  {end do\}
427  {end do\}
428  {end do\}
429 
430  end subroutine put_bc_comm_types
431 
432  !> do update ghost cells of all blocks including physical boundaries
433  subroutine getbc(time,qdt,psb,nwstart,nwbc,req_diag)
435  use mod_physics
436 
437  double precision, intent(in) :: time, qdt
438  type(state), target :: psb(max_blocks)
439  integer, intent(in) :: nwstart ! Fill from nwstart
440  integer, intent(in) :: nwbc ! Number of variables to fill
441  logical, intent(in), optional :: req_diag ! If false, skip diagonal ghost cells
442 
443  double precision :: time_bcin
444  integer :: ipole, nwhead, nwtail
445  integer :: iigrid, igrid, ineighbor, ipe_neighbor, isizes
446  integer :: ixR^L, ixS^L
447  integer :: i^D, n_i^D, ic^D, inc^D, n_inc^D, iib^D, idir
448  ! store physical boundary indicating index
449  integer :: idphyb(ndim,max_blocks)
450  integer :: isend_buf(npwbuf), ipwbuf, nghostcellsco
451  ! index pointer for buffer arrays as a start for a segment
452  integer :: ibuf_start, ibuf_next
453  ! shapes of reshape
454  integer, dimension(1) :: shapes
455  logical :: req_diagonal
456  type(wbuffer) :: pwbuf(npwbuf)
457 
458  time_bcin=mpi_wtime()
459 
460  nwhead=nwstart
461  nwtail=nwstart+nwbc-1
462 
463  req_diagonal = .true.
464  if (present(req_diag)) req_diagonal = req_diag
465 
466  ! fill internal physical boundary
467  if (internalboundary) then
468  call getintbc(time,ixg^ll)
469  end if
470 
471  ! fill physical-boundary ghost cells before internal ghost-cell values exchange
472  if(bcphys.and. .not.stagger_grid) then
473  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
474  do iigrid=1,igridstail; igrid=igrids(iigrid);
475  if(.not.phyboundblock(igrid)) cycle
476  call fill_boundary_before_gc(igrid)
477  end do
478  !$OMP END PARALLEL DO
479  end if
480 
481  ! prepare coarse values to send to coarser neighbors
482  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
483  do iigrid=1,igridstail; igrid=igrids(iigrid);
484  if(any(neighbor_type(:^d&,igrid)==neighbor_coarse)) then
485  call coarsen_grid(psb(igrid),ixg^ll,ixm^l,psc(igrid),ixcog^l,ixcom^l)
486  {do i^db=-1,1\}
487  if(skip_direction([ i^d ])) cycle
488  if(neighbor_type(i^d,igrid)==neighbor_coarse) call fill_coarse_boundary(igrid,i^d)
489  {end do\}
490  end if
491  end do
492  !$OMP END PARALLEL DO
493 
494  ! default : no singular axis
495  ipole=0
496  irecv_c=0
497  isend_c=0
498  isend_buf=0
499  ipwbuf=1
500 
501  if(stagger_grid) then
502  ibuf_recv_srl=1
503  ibuf_recv_r=1
504  ibuf_recv_p=1
505  ibuf_send_srl=1
506  ibuf_send_r=1
507  ibuf_send_p=1
508  irecv_srl=0
509  irecv_r=0
510  irecv_p=0
511  isend_srl=0
512  isend_r=0
513  isend_p=0
514  end if
515 
516  ! MPI receive ghost-cell values from sibling blocks and finer neighbors in different processors
517  do iigrid=1,igridstail; igrid=igrids(iigrid);
518  call identifyphysbound(ps(igrid),iib^d)
519  ^d&idphyb(^d,igrid)=iib^d;
520  {do i^db=-1,1\}
521  if (skip_direction([ i^d ])) cycle
522  select case (neighbor_type(i^d,igrid))
523  case (neighbor_sibling)
524  call bc_recv_srl
525  case (neighbor_fine)
526  call bc_recv_restrict
527  end select
528  {end do\}
529  end do
530 
531  ! MPI send ghost-cell values to sibling blocks and coarser neighbors in different processors
532  do iigrid=1,igridstail; igrid=igrids(iigrid);
533  ^d&iib^d=idphyb(^d,igrid);
534  {do i^db=-1,1\}
535  if(skip_direction([ i^d ])) cycle
536  select case (neighbor_type(i^d,igrid))
537  case (neighbor_sibling)
538  call bc_send_srl
539  case (neighbor_coarse)
540  call bc_send_restrict
541  end select
542  {end do\}
543  end do
544 
545  ! fill ghost-cell values of sibling blocks and coarser neighbors in the same processor
546  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid,iib^D)
547  do iigrid=1,igridstail; igrid=igrids(iigrid);
548  ^d&iib^d=idphyb(^d,igrid);
549  {do i^db=-1,1\}
550  if(skip_direction([ i^d ])) cycle
551  select case (neighbor_type(i^d,igrid))
552  case(neighbor_sibling)
553  call bc_fill_srl(igrid,i^d,iib^d)
554  case(neighbor_coarse)
555  call bc_fill_restrict(igrid,i^d,iib^d)
556  end select
557  {end do\}
558  end do
559  !$OMP END PARALLEL DO
560 
561  call mpi_waitall(irecv_c,recvrequest_c_sr,recvstatus_c_sr,ierrmpi)
562  call mpi_waitall(isend_c,sendrequest_c_sr,sendstatus_c_sr,ierrmpi)
563 
564  if(stagger_grid) then
565  call mpi_waitall(nrecv_bc_srl,recvrequest_srl,recvstatus_srl,ierrmpi)
566  call mpi_waitall(nsend_bc_srl,sendrequest_srl,sendstatus_srl,ierrmpi)
567  call mpi_waitall(nrecv_bc_r,recvrequest_r,recvstatus_r,ierrmpi)
568  call mpi_waitall(nsend_bc_r,sendrequest_r,sendstatus_r,ierrmpi)
569  ! unpack the received data from sibling blocks and finer neighbors to fill ghost-cell staggered values
570  ibuf_recv_srl=1
571  ibuf_recv_r=1
572  do iigrid=1,igridstail; igrid=igrids(iigrid);
573  ^d&iib^d=idphyb(^d,igrid);
574  {do i^db=-1,1\}
575  if (skip_direction([ i^d ])) cycle
576  select case (neighbor_type(i^d,igrid))
577  case (neighbor_sibling)
578  call bc_fill_srl_stg
579  case (neighbor_fine)
581  end select
582  {end do\}
583  end do
584  end if
585 
586  do ipwbuf=1,npwbuf
587  if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
588  end do
589 
590  irecv_c=0
591  isend_c=0
592  isend_buf=0
593  ipwbuf=1
594 
595  ! MPI receive ghost-cell values from coarser neighbors in different processors
596  do iigrid=1,igridstail; igrid=igrids(iigrid);
597  ^d&iib^d=idphyb(^d,igrid);
598  {do i^db=-1,1\}
599  if (skip_direction([ i^d ])) cycle
600  if (neighbor_type(i^d,igrid)==neighbor_coarse) call bc_recv_prolong
601  {end do\}
602  end do
603  ! MPI send ghost-cell values to finer neighbors in different processors
604  do iigrid=1,igridstail; igrid=igrids(iigrid);
605  ^d&iib^d=idphyb(^d,igrid);
606  {do i^db=-1,1\}
607  if (skip_direction([ i^d ])) cycle
608  if (neighbor_type(i^d,igrid)==neighbor_fine) call bc_send_prolong
609  {end do\}
610  end do
611 
612  ! fill coarse ghost-cell values of finer neighbors in the same processor
613  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid,iib^D)
614  do iigrid=1,igridstail; igrid=igrids(iigrid);
615  ^d&iib^d=idphyb(^d,igrid);
616  {do i^db=-1,1\}
617  if (skip_direction([ i^d ])) cycle
618  if (neighbor_type(i^d,igrid)==neighbor_fine) call bc_fill_prolong(igrid,i^d,iib^d)
619  {end do\}
620  end do
621  !$OMP END PARALLEL DO
622 
623  call mpi_waitall(irecv_c,recvrequest_c_p,recvstatus_c_p,ierrmpi)
624  call mpi_waitall(isend_c,sendrequest_c_p,sendstatus_c_p,ierrmpi)
625 
626  if(stagger_grid) then
627  call mpi_waitall(nrecv_bc_p,recvrequest_p,recvstatus_p,ierrmpi)
628  call mpi_waitall(nsend_bc_p,sendrequest_p,sendstatus_p,ierrmpi)
629 
630  ! fill coarser representative ghost cells after receipt
631  ibuf_recv_p=1
632  do iigrid=1,igridstail; igrid=igrids(iigrid);
633  ^d&iib^d=idphyb(^d,igrid);
634  {do i^db=-1,1\}
635  if (skip_direction([ i^d ])) cycle
636  if(neighbor_type(i^d,igrid)==neighbor_coarse) call bc_fill_prolong_stg
637  {end do\}
638  end do
639  end if
640  ! do prolongation on the ghost-cell values based on the received coarse values from coarser neighbors
641  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
642  do iigrid=1,igridstail; igrid=igrids(iigrid);
643  call gc_prolong(igrid)
644  end do
645  !$OMP END PARALLEL DO
646 
647  do ipwbuf=1,npwbuf
648  if (isend_buf(ipwbuf)/=0) deallocate(pwbuf(ipwbuf)%w)
649  end do
650 
651  ! fill physical boundary ghost cells after internal ghost-cell values exchange
652  if(bcphys.and.stagger_grid) then
653  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
654  do iigrid=1,igridstail; igrid=igrids(iigrid);
655  if(.not.phyboundblock(igrid)) cycle
656  call fill_boundary_after_gc(igrid)
657  end do
658  !$OMP END PARALLEL DO
659  end if
660 
661  ! modify normal component of magnetic field to fix divB=0
662  if(bcphys.and.associated(phys_boundary_adjust)) then
663  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
664  do iigrid=1,igridstail; igrid=igrids(iigrid);
665  if(.not.phyboundblock(igrid)) cycle
666  call phys_boundary_adjust(igrid,psb)
667  end do
668  !$OMP END PARALLEL DO
669  end if
670 
671  time_bc=time_bc+(mpi_wtime()-time_bcin)
672 
673  contains
674 
675  logical function skip_direction(dir)
676  integer, intent(in) :: dir(^nd)
677 
678  if (all(dir == 0)) then
679  skip_direction = .true.
680  else if (.not. req_diagonal .and. count(dir /= 0) > 1) then
681  skip_direction = .true.
682  else
683  skip_direction = .false.
684  end if
685  end function skip_direction
686 
687  !> Physical boundary conditions
688  subroutine fill_boundary_before_gc(igrid)
689 
690  integer, intent(in) :: igrid
691 
692  integer :: idims,iside,i^D,k^L,ixB^L
693 
694  block=>psb(igrid)
695  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
696  do idims=1,ndim
697  ! to avoid using as yet unknown corner info in more than 1D, we
698  ! fill only interior mesh ranges of the ghost cell ranges at first,
699  ! and progressively enlarge the ranges to include corners later
700  {
701  kmin^d=merge(0, 1, idims==^d)
702  kmax^d=merge(0, 1, idims==^d)
703  ixbmin^d=ixglo^d+kmin^d*nghostcells
704  ixbmax^d=ixghi^d-kmax^d*nghostcells
705  \}
706  {^iftwod
707  if(idims > 1 .and. neighbor_type(-1,0,igrid)==neighbor_boundary) ixbmin1=ixglo1
708  if(idims > 1 .and. neighbor_type( 1,0,igrid)==neighbor_boundary) ixbmax1=ixghi1}
709  {^ifthreed
710  if(idims > 1 .and. neighbor_type(-1,0,0,igrid)==neighbor_boundary) ixbmin1=ixglo1
711  if(idims > 1 .and. neighbor_type( 1,0,0,igrid)==neighbor_boundary) ixbmax1=ixghi1
712  if(idims > 2 .and. neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixbmin2=ixglo2
713  if(idims > 2 .and. neighbor_type(0, 1,0,igrid)==neighbor_boundary) ixbmax2=ixghi2}
714  do iside=1,2
715  i^d=kr(^d,idims)*(2*iside-3);
716  if (aperiodb(idims)) then
717  if (neighbor_type(i^d,igrid) /= neighbor_boundary .and. &
718  .not. psb(igrid)%is_physical_boundary(2*idims-2+iside)) cycle
719  else
720  if (neighbor_type(i^d,igrid) /= neighbor_boundary) cycle
721  end if
722  call bc_phys(iside,idims,time,qdt,psb(igrid),ixg^ll,ixb^l)
723  end do
724  end do
725 
726  end subroutine fill_boundary_before_gc
727 
728  !> Physical boundary conditions
729  subroutine fill_boundary_after_gc(igrid)
730 
731  integer, intent(in) :: igrid
732 
733  integer :: idims,iside,i^D,k^L,ixB^L
734 
735  block=>psb(igrid)
736  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
737  do idims=1,ndim
738  ! to avoid using as yet unknown corner info in more than 1D, we
739  ! fill only interior mesh ranges of the ghost cell ranges at first,
740  ! and progressively enlarge the ranges to include corners later
741  kmin1=0; kmax1=0;
742  {^iftwod
743  kmin2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0,-1,igrid)==1)
744  kmax2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0, 1,igrid)==1)}
745  {^ifthreed
746  kmin2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0,-1,0,igrid)==1)
747  kmax2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0, 1,0,igrid)==1)
748  kmin3=merge(1, 0, idims .lt. 3 .and. neighbor_type(0,0,-1,igrid)==1)
749  kmax3=merge(1, 0, idims .lt. 3 .and. neighbor_type(0,0, 1,igrid)==1)}
750  ixbmin^d=ixglo^d+kmin^d*nghostcells;
751  ixbmax^d=ixghi^d-kmax^d*nghostcells;
752  do iside=1,2
753  i^d=kr(^d,idims)*(2*iside-3);
754  if (aperiodb(idims)) then
755  if (neighbor_type(i^d,igrid) /= neighbor_boundary .and. &
756  .not. psb(igrid)%is_physical_boundary(2*idims-2+iside)) cycle
757  else
758  if (neighbor_type(i^d,igrid) /= neighbor_boundary) cycle
759  end if
760  call bc_phys(iside,idims,time,qdt,psb(igrid),ixg^ll,ixb^l)
761  end do
762  end do
763 
764  end subroutine fill_boundary_after_gc
765 
766  !> Receive from sibling at same refinement level
767  subroutine bc_recv_srl
768 
769  ipe_neighbor=neighbor(2,i^d,igrid)
770  if (ipe_neighbor/=mype) then
771  irecv_c=irecv_c+1
772  itag=(3**^nd+4**^nd)*(igrid-1)+{(i^d+1)*3**(^d-1)+}
773  call mpi_irecv(psb(igrid)%w,1,type_recv_srl(iib^d,i^d), &
774  ipe_neighbor,itag,icomm,recvrequest_c_sr(irecv_c),ierrmpi)
775  if(stagger_grid) then
777  call mpi_irecv(recvbuffer_srl(ibuf_recv_srl),sizes_srl_recv_total(i^d),mpi_double_precision, &
778  ipe_neighbor,itag,icomm,recvrequest_srl(irecv_srl),ierrmpi)
780  end if
781  end if
782 
783  end subroutine bc_recv_srl
784 
785  !> Receive from fine neighbor
786  subroutine bc_recv_restrict
787 
788  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
789  inc^db=2*i^db+ic^db\}
790  ipe_neighbor=neighbor_child(2,inc^d,igrid)
791  if (ipe_neighbor/=mype) then
792  irecv_c=irecv_c+1
793  itag=(3**^nd+4**^nd)*(igrid-1)+3**^nd+{inc^d*4**(^d-1)+}
794  call mpi_irecv(psb(igrid)%w,1,type_recv_r(iib^d,inc^d), &
795  ipe_neighbor,itag,icomm,recvrequest_c_sr(irecv_c),ierrmpi)
796  if(stagger_grid) then
797  irecv_r=irecv_r+1
798  call mpi_irecv(recvbuffer_r(ibuf_recv_r),sizes_r_recv_total(inc^d), &
799  mpi_double_precision,ipe_neighbor,itag, &
800  icomm,recvrequest_r(irecv_r),ierrmpi)
802  end if
803  end if
804  {end do\}
805 
806  end subroutine bc_recv_restrict
807 
808  !> Send to sibling at same refinement level
809  subroutine bc_send_srl
810 
811  ipe_neighbor=neighbor(2,i^d,igrid)
812 
813  if(ipe_neighbor/=mype) then
814  ineighbor=neighbor(1,i^d,igrid)
815  ipole=neighbor_pole(i^d,igrid)
816  if(ipole==0) then
817  n_i^d=-i^d;
818  isend_c=isend_c+1
819  itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
820  call mpi_isend(psb(igrid)%w,1,type_send_srl(iib^d,i^d), &
821  ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
822  if(stagger_grid) then
823  ibuf_start=ibuf_send_srl
824  do idir=1,ndim
825  ixs^l=ixs_srl_stg_^l(idir,i^d);
826  ibuf_next=ibuf_start+sizes_srl_send_stg(idir,i^d)
827  shapes=(/sizes_srl_send_stg(idir,i^d)/)
828  sendbuffer_srl(ibuf_start:ibuf_next-1)=&
829  reshape(psb(igrid)%ws(ixs^s,idir),shapes)
830  ibuf_start=ibuf_next
831  end do
833  call mpi_isend(sendbuffer_srl(ibuf_send_srl),sizes_srl_send_total(i^d),&
834  mpi_double_precision, ipe_neighbor,itag,icomm,sendrequest_srl(isend_srl),ierrmpi)
835  ibuf_send_srl=ibuf_next
836  end if
837  else
838  ixs^l=ixs_srl_^l(iib^d,i^d);
839  select case (ipole)
840  {case (^d)
841  n_i^d=i^d^d%n_i^dd=-i^dd;\}
842  end select
843  if (isend_buf(ipwbuf)/=0) then
844  call mpi_wait(sendrequest_c_sr(isend_buf(ipwbuf)), &
845  sendstatus_c_sr(:,isend_buf(ipwbuf)),ierrmpi)
846  deallocate(pwbuf(ipwbuf)%w)
847  end if
848  allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
849  call pole_buffer(pwbuf(ipwbuf)%w,ixs^l,ixs^l,psb(igrid)%w,ixg^ll,ixs^l)
850  isend_c=isend_c+1
851  isend_buf(ipwbuf)=isend_c
852  itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
853  isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
854  call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
855  ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
856  ipwbuf=1+modulo(ipwbuf,npwbuf)
857  if(stagger_grid) then
858  ibuf_start=ibuf_send_srl
859  do idir=1,ndim
860  ixs^l=ixs_srl_stg_^l(idir,i^d);
861  ibuf_next=ibuf_start+sizes_srl_send_stg(idir,i^d)
862  shapes=(/sizes_srl_send_stg(idir,i^d)/)
863  sendbuffer_srl(ibuf_start:ibuf_next-1)=&
864  reshape(psb(igrid)%ws(ixs^s,idir),shapes)
865  ibuf_start=ibuf_next
866  end do
868  call mpi_isend(sendbuffer_srl(ibuf_send_srl),sizes_srl_send_total(i^d),&
869  mpi_double_precision, ipe_neighbor,itag,icomm,sendrequest_srl(isend_srl),ierrmpi)
870  ibuf_send_srl=ibuf_next
871  end if
872  end if
873  end if
874 
875  end subroutine bc_send_srl
876 
877  subroutine bc_fill_srl(igrid,i^D,iib^D)
878  integer, intent(in) :: igrid,i^D,iib^D
879  integer :: ineighbor,ipe_neighbor,ipole,ixS^L,ixR^L,n_i^D,idir
880 
881  ipe_neighbor=neighbor(2,i^d,igrid)
882  if(ipe_neighbor==mype) then
883  ineighbor=neighbor(1,i^d,igrid)
884  ipole=neighbor_pole(i^d,igrid)
885  if(ipole==0) then
886  n_i^d=-i^d;
887  ixs^l=ixs_srl_^l(iib^d,i^d);
888  ixr^l=ixr_srl_^l(iib^d,n_i^d);
889  psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
890  psb(igrid)%w(ixs^s,nwhead:nwtail)
891  if(stagger_grid) then
892  do idir=1,ndim
893  ixs^l=ixs_srl_stg_^l(idir,i^d);
894  ixr^l=ixr_srl_stg_^l(idir,n_i^d);
895  psb(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
896  end do
897  end if
898  else
899  ixs^l=ixs_srl_^l(iib^d,i^d);
900  select case (ipole)
901  {case (^d)
902  n_i^d=i^d^d%n_i^dd=-i^dd;\}
903  end select
904  ixr^l=ixr_srl_^l(iib^d,n_i^d);
905  call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^l,psb(igrid)%w,ixg^ll,ixs^l,ipole)
906  if(stagger_grid) then
907  do idir=1,ndim
908  ixs^l=ixs_srl_stg_^l(idir,i^d);
909  ixr^l=ixr_srl_stg_^l(idir,n_i^d);
910  call pole_copy_stg(psb(ineighbor)%ws,ixgs^ll,ixr^l,psb(igrid)%ws,ixgs^ll,ixs^l,idir,ipole)
911  end do
912  end if
913  end if
914  end if
915 
916  end subroutine bc_fill_srl
917 
918  subroutine fill_coarse_boundary(igrid,i^D)
919  integer, intent(in) :: igrid,i^D
920 
921  integer :: idims,iside,k^L,ixB^L,ii^D
922 
923  if(phyboundblock(igrid).and..not.stagger_grid.and.bcphys) then
924  ! to use block in physical boundary setup for coarse representative
925  block=>psc(igrid)
926  ! filling physical boundary ghost cells of a coarser representative block for
927  ! sending swap region with width of nghostcells to its coarser neighbor
928  do idims=1,ndim
929  ! to avoid using as yet unknown corner info in more than 1D, we
930  ! fill only interior mesh ranges of the ghost cell ranges at first,
931  ! and progressively enlarge the ranges to include corners later
932  {kmin^d=merge(0, 1, idims==^d)
933  kmax^d=merge(0, 1, idims==^d)
934  ixbmin^d=ixcogmin^d+kmin^d*nghostcells
935  ixbmax^d=ixcogmax^d-kmax^d*nghostcells\}
936  {^iftwod
937  if(idims > 1 .and. neighbor_type(-1,0,igrid)==neighbor_boundary) ixbmin1=ixcogmin1
938  if(idims > 1 .and. neighbor_type( 1,0,igrid)==neighbor_boundary) ixbmax1=ixcogmax1}
939  {^ifthreed
940  if(idims > 1 .and. neighbor_type(-1,0,0,igrid)==neighbor_boundary) ixbmin1=ixcogmin1
941  if(idims > 1 .and. neighbor_type( 1,0,0,igrid)==neighbor_boundary) ixbmax1=ixcogmax1
942  if(idims > 2 .and. neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixbmin2=ixcogmin2
943  if(idims > 2 .and. neighbor_type(0, 1,0,igrid)==neighbor_boundary) ixbmax2=ixcogmax2}
944  {if(i^d==-1) then
945  ixbmin^d=ixcogmin^d+nghostcells
946  ixbmax^d=ixcogmin^d+2*nghostcells-1
947  else if(i^d==1) then
948  ixbmin^d=ixcogmax^d-2*nghostcells+1
949  ixbmax^d=ixcogmax^d-nghostcells
950  end if\}
951  do iside=1,2
952  ii^d=kr(^d,idims)*(2*iside-3);
953  if ({abs(i^d)==1.and.abs(ii^d)==1|.or.}) cycle
954  if (neighbor_type(ii^d,igrid)/=neighbor_boundary) cycle
955  call bc_phys(iside,idims,time,0.d0,psc(igrid),ixcog^l,ixb^l)
956  end do
957  end do
958  end if
959 
960  end subroutine fill_coarse_boundary
961 
962  !> Send to coarser neighbor
963  subroutine bc_send_restrict
964 
965  ipe_neighbor=neighbor(2,i^d,igrid)
966  if(ipe_neighbor/=mype) then
967  ic^d=1+modulo(node(pig^d_,igrid)-1,2);
968  if({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
969  ineighbor=neighbor(1,i^d,igrid)
970  ipole=neighbor_pole(i^d,igrid)
971  if(ipole==0) then
972  n_inc^d=-2*i^d+ic^d;
973  isend_c=isend_c+1
974  itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
975  call mpi_isend(psc(igrid)%w,1,type_send_r(iib^d,i^d), &
976  ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
977  if(stagger_grid) then
978  ibuf_start=ibuf_send_r
979  do idir=1,ndim
980  ixs^l=ixs_r_stg_^l(idir,i^d);
981  ibuf_next=ibuf_start+sizes_r_send_stg(idir,i^d)
982  shapes=(/sizes_r_send_stg(idir,i^d)/)
983  sendbuffer_r(ibuf_start:ibuf_next-1)=&
984  reshape(psc(igrid)%ws(ixs^s,idir),shapes)
985  ibuf_start=ibuf_next
986  end do
987  isend_r=isend_r+1
988  call mpi_isend(sendbuffer_r(ibuf_send_r),sizes_r_send_total(i^d),&
989  mpi_double_precision,ipe_neighbor,itag, &
990  icomm,sendrequest_r(isend_r),ierrmpi)
991  ibuf_send_r=ibuf_next
992  end if
993  else
994  ixs^l=ixs_r_^l(iib^d,i^d);
995  select case (ipole)
996  {case (^d)
997  n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
998  end select
999  if(isend_buf(ipwbuf)/=0) then
1000  call mpi_wait(sendrequest_c_sr(isend_buf(ipwbuf)), &
1001  sendstatus_c_sr(:,isend_buf(ipwbuf)),ierrmpi)
1002  deallocate(pwbuf(ipwbuf)%w)
1003  end if
1004  allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
1005  call pole_buffer(pwbuf(ipwbuf)%w,ixs^l,ixs^l,psc(igrid)%w,ixcog^l,ixs^l)
1006  isend_c=isend_c+1
1007  isend_buf(ipwbuf)=isend_c
1008  itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
1009  isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
1010  call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
1011  ipe_neighbor,itag,icomm,sendrequest_c_sr(isend_c),ierrmpi)
1012  ipwbuf=1+modulo(ipwbuf,npwbuf)
1013  if(stagger_grid) then
1014  ibuf_start=ibuf_send_r
1015  do idir=1,ndim
1016  ixs^l=ixs_r_stg_^l(idir,i^d);
1017  ibuf_next=ibuf_start+sizes_r_send_stg(idir,i^d)
1018  shapes=(/sizes_r_send_stg(idir,i^d)/)
1019  sendbuffer_r(ibuf_start:ibuf_next-1)=&
1020  reshape(psc(igrid)%ws(ixs^s,idir),shapes)
1021  ibuf_start=ibuf_next
1022  end do
1023  isend_r=isend_r+1
1024  call mpi_isend(sendbuffer_r(ibuf_send_r),sizes_r_send_total(i^d),&
1025  mpi_double_precision,ipe_neighbor,itag, &
1026  icomm,sendrequest_r(isend_r),ierrmpi)
1027  ibuf_send_r=ibuf_next
1028  end if
1029  end if
1030  end if
1031 
1032  end subroutine bc_send_restrict
1033 
1034  !> fill coarser neighbor's ghost cells
1035  subroutine bc_fill_restrict(igrid,i^D,iib^D)
1036  integer, intent(in) :: igrid,i^D,iib^D
1037 
1038  integer :: ic^D,n_inc^D,ixS^L,ixR^L,ipe_neighbor,ineighbor,ipole,idir
1039 
1040  ipe_neighbor=neighbor(2,i^d,igrid)
1041  if(ipe_neighbor==mype) then
1042  ic^d=1+modulo(node(pig^d_,igrid)-1,2);
1043  if({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
1044  ineighbor=neighbor(1,i^d,igrid)
1045  ipole=neighbor_pole(i^d,igrid)
1046  if(ipole==0) then
1047  n_inc^d=-2*i^d+ic^d;
1048  ixs^l=ixs_r_^l(iib^d,i^d);
1049  ixr^l=ixr_r_^l(iib^d,n_inc^d);
1050  psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
1051  psc(igrid)%w(ixs^s,nwhead:nwtail)
1052  if(stagger_grid) then
1053  do idir=1,ndim
1054  ixs^l=ixs_r_stg_^l(idir,i^d);
1055  ixr^l=ixr_r_stg_^l(idir,n_inc^d);
1056  psb(ineighbor)%ws(ixr^s,idir)=psc(igrid)%ws(ixs^s,idir)
1057  end do
1058  end if
1059  else
1060  ixs^l=ixs_r_^l(iib^d,i^d);
1061  select case (ipole)
1062  {case (^d)
1063  n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
1064  end select
1065  ixr^l=ixr_r_^l(iib^d,n_inc^d);
1066  call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^l,psc(igrid)%w,ixcog^l,ixs^l,ipole)
1067  if(stagger_grid) then
1068  do idir=1,ndim
1069  ixs^l=ixs_r_stg_^l(idir,i^d);
1070  ixr^l=ixr_r_stg_^l(idir,n_inc^d);
1071  !! Fill ghost cells
1072  call pole_copy_stg(psb(ineighbor)%ws,ixgs^ll,ixr^l,psc(igrid)%ws,ixcogs^l,ixs^l,idir,ipole)
1073  end do
1074  end if
1075  end if
1076  end if
1077 
1078  end subroutine bc_fill_restrict
1079 
1080  !> fill siblings ghost cells with received data
1081  subroutine bc_fill_srl_stg
1082  double precision :: tmp(ixGs^T)
1083  integer :: ixS^L,ixR^L,n_i^D,ixSsync^L,ixRsync^L
1084  integer :: idir
1085 
1086  ipe_neighbor=neighbor(2,i^d,igrid)
1087  if(ipe_neighbor/=mype) then
1088  ineighbor=neighbor(1,i^d,igrid)
1089  ipole=neighbor_pole(i^d,igrid)
1090 
1091  !! Now the special treatment of the pole is done here, at the receive step
1092  if (ipole==0) then
1093  ixr^l=ixr_srl_^l(iib^d,i^d);
1094  !! Unpack the buffer and fill the ghost cells
1095  n_i^d=-i^d;
1096  do idir=1,ndim
1097  ixs^l=ixs_srl_stg_^l(idir,n_i^d);
1098  ixr^l=ixr_srl_stg_^l(idir,i^d);
1099  ibuf_next=ibuf_recv_srl+sizes_srl_recv_stg(idir,i^d)
1100  tmp(ixs^s) = reshape(source=recvbuffer_srl(ibuf_recv_srl:ibuf_next-1),shape=shape(psb(igrid)%ws(ixs^s,idir)))
1101  psb(igrid)%ws(ixr^s,idir) = tmp(ixs^s)
1102  ibuf_recv_srl=ibuf_next
1103  end do
1104  else ! There is a pole
1105  select case (ipole)
1106  {case (^d)
1107  n_i^d=i^d^d%n_i^dd=-i^dd;\}
1108  end select
1109  pole_buf%ws=zero
1110  do idir=1,ndim
1111  ixr^l=ixr_srl_stg_^l(idir,i^d);
1112  ixs^l=ixs_srl_stg_^l(idir,n_i^d);
1113  ibuf_next=ibuf_recv_srl+sizes_srl_recv_stg(idir,i^d)
1114  pole_buf%ws(ixs^s,idir)=reshape(source=recvbuffer_srl(ibuf_recv_srl:ibuf_next-1),&
1115  shape=shape(psb(igrid)%ws(ixs^s,idir)))
1116  ibuf_recv_srl=ibuf_next
1117  call pole_copy_stg(psb(igrid)%ws,ixgs^ll,ixr^l,pole_buf%ws,ixgs^ll,ixs^l,idir,ipole)
1118  end do
1119  end if
1120  end if
1121 
1122  end subroutine bc_fill_srl_stg
1123 
1124  subroutine indices_for_syncing(idir,i^D,ixR^L,ixS^L,ixRsync^L,ixSsync^L)
1125  integer, intent(in) :: i^D,idir
1126  integer, intent(inout) :: ixR^L,ixS^L
1127  integer, intent(out) :: ixRsync^L,ixSsync^L
1128 
1129  ixrsync^l=ixr^l;
1130  ixssync^l=ixs^l;
1131 
1132  {
1133  if (i^d == -1 .and. idir == ^d) then
1134  ixrsyncmin^d = ixrmax^d
1135  ixrsyncmax^d = ixrmax^d
1136  ixssyncmin^d = ixsmax^d
1137  ixssyncmax^d = ixsmax^d
1138  ixrmax^d = ixrmax^d - 1
1139  ixsmax^d = ixsmax^d - 1
1140  else if (i^d == 1 .and. idir == ^d) then
1141  ixrsyncmin^d = ixrmin^d
1142  ixrsyncmax^d = ixrmin^d
1143  ixssyncmin^d = ixsmin^d
1144  ixssyncmax^d = ixsmin^d
1145  ixrmin^d = ixrmin^d + 1
1146  ixsmin^d = ixsmin^d + 1
1147  end if
1148  \}
1149 
1150  end subroutine indices_for_syncing
1151 
1152  !> fill restricted ghost cells after receipt
1154 
1155  ipole=neighbor_pole(i^d,igrid)
1156  if (ipole==0) then
1157  ! Loop over the children ic^D to and their neighbors inc^D
1158  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1159  inc^db=2*i^db+ic^db\}
1160  ipe_neighbor=neighbor_child(2,inc^d,igrid)
1161  if(ipe_neighbor/=mype) then
1162  ineighbor=neighbor_child(1,inc^d,igrid)
1163  n_i^d=-i^d;
1164  !! Unpack the buffer and fill the ghost cells
1165  do idir=1,ndim
1166  ixr^l=ixr_r_stg_^l(idir,inc^d);
1167  ibuf_next=ibuf_recv_r+sizes_r_recv_stg(idir,inc^d)
1168  psb(igrid)%ws(ixr^s,idir)=reshape(source=recvbuffer_r(ibuf_recv_r:ibuf_next-1),&
1169  shape=shape(psb(igrid)%ws(ixr^s,idir)))
1170  ibuf_recv_r=ibuf_next
1171  end do
1172  end if
1173  {end do\}
1174  else !! There is a pole
1175  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1176  inc^db=2*i^db+ic^db\}
1177  ipe_neighbor=neighbor_child(2,inc^d,igrid)
1178  if(ipe_neighbor/=mype) then
1179  ineighbor=neighbor_child(1,inc^d,igrid)
1180  select case(ipole)
1181  {case (^d)
1182  n_i^d=i^d^d%n_i^dd=-i^dd;\}
1183  end select
1184  ixr^l=ixr_r_^l(iib^d,inc^d);
1185  !! Unpack the buffer and fill an auxiliary array
1186  pole_buf%ws=zero
1187  do idir=1,ndim
1188  ixs^l=ixs_r_stg_^l(idir,n_i^d);
1189  ixr^l=ixr_r_stg_^l(idir,inc^d);
1190  ibuf_next=ibuf_recv_r+sizes_r_recv_stg(idir,inc^d)
1191  pole_buf%ws(ixr^s,idir)=reshape(source=recvbuffer_r(ibuf_recv_r:ibuf_next-1),&
1192  shape=shape(psb(igrid)%ws(ixr^s,idir)))
1193  call pole_copy_stg(psb(igrid)%ws,ixgs^ll,ixr^l,pole_buf%ws,ixgs^ll,ixr^l,idir,ipole)
1194  ibuf_recv_r=ibuf_next
1195  end do
1196  end if
1197  {end do\}
1198  end if
1199 
1200  end subroutine bc_fill_restrict_stg
1201 
1202  !> Receive from coarse neighbor
1203  subroutine bc_recv_prolong
1204 
1205  ic^d=1+modulo(node(pig^d_,igrid)-1,2);
1206  if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
1207 
1208  ipe_neighbor=neighbor(2,i^d,igrid)
1209  if (ipe_neighbor/=mype) then
1210  irecv_c=irecv_c+1
1211  inc^d=ic^d+i^d;
1212  itag=(3**^nd+4**^nd)*(igrid-1)+3**^nd+{inc^d*4**(^d-1)+}
1213  call mpi_irecv(psc(igrid)%w,1,type_recv_p(iib^d,inc^d), &
1214  ipe_neighbor,itag,icomm,recvrequest_c_p(irecv_c),ierrmpi)
1215  if(stagger_grid) then
1216  irecv_p=irecv_p+1
1217  call mpi_irecv(recvbuffer_p(ibuf_recv_p),sizes_p_recv_total(inc^d),&
1218  mpi_double_precision,ipe_neighbor,itag,&
1219  icomm,recvrequest_p(irecv_p),ierrmpi)
1221  end if
1222  end if
1223 
1224  end subroutine bc_recv_prolong
1225 
1226  !> Send to finer neighbor
1227  subroutine bc_send_prolong
1228  integer :: ii^D
1229 
1230  ipole=neighbor_pole(i^d,igrid)
1231 
1232  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1233  inc^db=2*i^db+ic^db\}
1234  ipe_neighbor=neighbor_child(2,inc^d,igrid)
1235  if(ipe_neighbor/=mype) then
1236  ixs^l=ixs_p_^l(iib^d,inc^d);
1237  ineighbor=neighbor_child(1,inc^d,igrid)
1238  if(ipole==0) then
1239  n_i^d=-i^d;
1240  n_inc^d=ic^d+n_i^d;
1241  isend_c=isend_c+1
1242  itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
1243  call mpi_isend(psb(igrid)%w,1,type_send_p(iib^d,inc^d), &
1244  ipe_neighbor,itag,icomm,sendrequest_c_p(isend_c),ierrmpi)
1245  if(stagger_grid) then
1246  ibuf_start=ibuf_send_p
1247  do idir=1,ndim
1248  ixs^l=ixs_p_stg_^l(idir,inc^d);
1249  ibuf_next=ibuf_start+sizes_p_send_stg(idir,inc^d)
1250  shapes=(/sizes_p_send_stg(idir,inc^d)/)
1251  sendbuffer_p(ibuf_start:ibuf_next-1)=&
1252  reshape(psb(igrid)%ws(ixs^s,idir),shapes)
1253  ibuf_start=ibuf_next
1254  end do
1255  isend_p=isend_p+1
1256  call mpi_isend(sendbuffer_p(ibuf_send_p),sizes_p_send_total(inc^d),&
1257  mpi_double_precision,ipe_neighbor,itag, &
1258  icomm,sendrequest_p(isend_p),ierrmpi)
1259  ibuf_send_p=ibuf_next
1260  end if
1261  else
1262  select case (ipole)
1263  {case (^d)
1264  n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
1265  end select
1266  if(isend_buf(ipwbuf)/=0) then
1267  call mpi_wait(sendrequest_c_p(isend_buf(ipwbuf)), &
1268  sendstatus_c_p(:,isend_buf(ipwbuf)),ierrmpi)
1269  deallocate(pwbuf(ipwbuf)%w)
1270  end if
1271  allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
1272  call pole_buffer(pwbuf(ipwbuf)%w,ixs^l,ixs^l,psb(igrid)%w,ixg^ll,ixs^l)
1273  isend_c=isend_c+1
1274  isend_buf(ipwbuf)=isend_c
1275  itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
1276  isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
1277  call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
1278  ipe_neighbor,itag,icomm,sendrequest_c_p(isend_c),ierrmpi)
1279  ipwbuf=1+modulo(ipwbuf,npwbuf)
1280  if(stagger_grid) then
1281  ibuf_start=ibuf_send_p
1282  do idir=1,ndim
1283  ixs^l=ixs_p_stg_^l(idir,inc^d);
1284  ibuf_next=ibuf_start+sizes_p_send_stg(idir,inc^d)
1285  shapes=(/sizes_p_send_stg(idir,inc^d)/)
1286  sendbuffer_p(ibuf_start:ibuf_next-1)=&
1287  reshape(psb(igrid)%ws(ixs^s,idir),shapes)
1288  ibuf_start=ibuf_next
1289  end do
1290  isend_p=isend_p+1
1291  call mpi_isend(sendbuffer_p(ibuf_send_p),sizes_p_send_total(inc^d),&
1292  mpi_double_precision,ipe_neighbor,itag, &
1293  icomm,sendrequest_p(isend_p),ierrmpi)
1294  ibuf_send_p=ibuf_next
1295  end if
1296  end if
1297  end if
1298  {end do\}
1299 
1300  end subroutine bc_send_prolong
1301 
1302  !> Send to finer neighbor
1303  subroutine bc_fill_prolong(igrid,i^D,iib^D)
1304  integer, intent(in) :: igrid,i^D,iib^D
1305 
1306  integer :: ipe_neighbor,ineighbor,ixS^L,ixR^L,ic^D,inc^D,ipole,idir
1307 
1308  ipole=neighbor_pole(i^d,igrid)
1309 
1310  if(ipole==0) then
1311  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1312  inc^db=2*i^db+ic^db\}
1313  ipe_neighbor=neighbor_child(2,inc^d,igrid)
1314  if(ipe_neighbor==mype) then
1315  ixs^l=ixs_p_^l(iib^d,inc^d);
1316  ineighbor=neighbor_child(1,inc^d,igrid)
1317  ipole=neighbor_pole(i^d,igrid)
1318  n_i^d=-i^d;
1319  n_inc^d=ic^d+n_i^d;
1320  ixr^l=ixr_p_^l(iib^d,n_inc^d);
1321  psc(ineighbor)%w(ixr^s,nwhead:nwtail) &
1322  =psb(igrid)%w(ixs^s,nwhead:nwtail)
1323  if(stagger_grid) then
1324  do idir=1,ndim
1325  ixs^l=ixs_p_stg_^l(idir,inc^d);
1326  ixr^l=ixr_p_stg_^l(idir,n_inc^d);
1327  psc(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
1328  end do
1329  end if
1330  end if
1331  {end do\}
1332  else
1333  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1334  inc^db=2*i^db+ic^db\}
1335  ipe_neighbor=neighbor_child(2,inc^d,igrid)
1336  if(ipe_neighbor==mype) then
1337  ixs^l=ixs_p_^l(iib^d,inc^d);
1338  ineighbor=neighbor_child(1,inc^d,igrid)
1339  ipole=neighbor_pole(i^d,igrid)
1340  select case (ipole)
1341  {case (^d)
1342  n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
1343  end select
1344  ixr^l=ixr_p_^l(iib^d,n_inc^d);
1345  call pole_copy(psc(ineighbor)%w,ixcog^l,ixr^l,psb(igrid)%w,ixg^ll,ixs^l,ipole)
1346  if(stagger_grid) then
1347  do idir=1,ndim
1348  ixs^l=ixs_p_stg_^l(idir,inc^d);
1349  ixr^l=ixr_p_stg_^l(idir,n_inc^d);
1350  call pole_copy_stg(psc(ineighbor)%ws,ixcogs^l,ixr^l,psb(igrid)%ws,ixgs^ll,ixs^l,idir,ipole)
1351  end do
1352  end if
1353  end if
1354  {end do\}
1355  end if
1356  end subroutine bc_fill_prolong
1357 
1358  subroutine gc_prolong(igrid)
1359  integer, intent(in) :: igrid
1360 
1361  integer :: iib^D,i^D,idims,iside
1362  logical,dimension(-1:1^D&) :: NeedProlong
1363 
1364  ^d&iib^d=idphyb(^d,igrid);
1365  needprolong=.false.
1366  {do i^db=-1,1\}
1367  if (skip_direction([ i^d ])) cycle
1368  if (neighbor_type(i^d,igrid)==neighbor_coarse) then
1369  call bc_prolong(igrid,i^d,iib^d)
1370  needprolong(i^d)=.true.
1371  end if
1372  {end do\}
1373  if(stagger_grid) then
1374  ! Ghost cell prolongation for staggered variables
1375  ! must be done in a specific order.
1376  ! First the first neighbours, which have 2 indices=0 in 3D
1377  ! or one index=0 in 2D
1378  block=>psb(igrid)
1379  do idims=1,ndim
1380  i^d=0;
1381  select case(idims)
1382  {case(^d)
1383  do i^d=-1,1,2
1384  if (needprolong(i^dd)) call bc_prolong_stg(igrid,i^dd,iib^dd,needprolong)
1385  end do
1386  \}
1387  end select
1388  end do
1389  ! Then the second neighbours which have 1 index=0 in 3D
1390  ! (Only in 3D)
1391  {^ifthreed
1392  i1=0;
1393  do i2=-1,1,2
1394  do i3=-1,1,2
1395  if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
1396  end do
1397  end do
1398  i2=0;
1399  do i3=-1,1,2
1400  do i1=-1,1,2
1401  if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
1402  end do
1403  end do
1404  i3=0;
1405  do i1=-1,1,2
1406  do i2=-1,1,2
1407  if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
1408  end do
1409  end do
1410  }
1411  ! Finally, the corners, that have no index=0
1412  {do i^d=-1,1,2\}
1413  if (needprolong(i^d)) call bc_prolong_stg(igrid,i^d,iib^d,needprolong)
1414  {end do\}
1415  end if
1416  end subroutine gc_prolong
1417 
1418  !> fill coarser representative with data from coarser neighbors
1420  ic^d=1+modulo(node(pig^d_,igrid)-1,2);
1421  if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.}) return
1422 
1423  ipe_neighbor=neighbor(2,i^d,igrid)
1424  if(ipe_neighbor/=mype) then
1425  ineighbor=neighbor(1,i^d,igrid)
1426  ipole=neighbor_pole(i^d,igrid)
1427 
1428  if (ipole==0) then !! There is no pole
1429  inc^d=ic^d+i^d;
1430  ixr^l=ixr_p_^l(iib^d,inc^d);
1431  do idir=1,ndim
1432  ixr^l=ixr_p_stg_^l(idir,inc^d);
1433  ibuf_next=ibuf_recv_p+sizes_p_recv_stg(idir,inc^d)
1434  psc(igrid)%ws(ixr^s,idir)=reshape(source=recvbuffer_p(ibuf_recv_p:ibuf_next-1),&
1435  shape=shape(psc(igrid)%ws(ixr^s,idir)))
1436  ibuf_recv_p=ibuf_next
1437  end do
1438  else !! There is a pole
1439  inc^d=ic^d+i^d;
1440  select case (ipole)
1441  {case (^d)
1442  n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
1443  end select
1444  !! Unpack the buffer and fill an auxiliary array
1445  pole_buf%ws=zero
1446  do idir=1,ndim
1447  ixr^l=ixr_p_stg_^l(idir,inc^d);
1448  ibuf_next=ibuf_recv_p+sizes_p_recv_stg(idir,inc^d)
1449  pole_buf%ws(ixr^s,idir)=reshape(source=recvbuffer_p(ibuf_recv_p:ibuf_next-1),&
1450  shape=shape(psc(igrid)%ws(ixr^s,idir)))
1451  call pole_copy_stg(psc(igrid)%ws,ixcogs^l,ixr^l,pole_buf%ws,ixgs^ll,ixr^l,idir,ipole)
1452  ibuf_recv_p=ibuf_next
1453  end do
1454  end if
1455  end if
1456 
1457  end subroutine bc_fill_prolong_stg
1458 
1459  !> do prolongation for fine blocks after receipt data from coarse neighbors
1460  subroutine bc_prolong(igrid,i^D,iib^D)
1462 
1463  integer :: i^D,iib^D,igrid
1464  integer :: ixFi^L,ixCo^L,ii^D, idims,iside,ixB^L
1465  double precision :: dxFi^D, dxCo^D, xFimin^D, xComin^D, invdxCo^D
1466 
1467  ixfi^l=ixr_srl_^l(iib^d,i^d);
1468  dxfi^d=rnode(rpdx^d_,igrid);
1469  dxco^d=two*dxfi^d;
1470  invdxco^d=1.d0/dxco^d;
1471 
1472  ! compute the enlarged grid lower left corner coordinates
1473  ! these are true coordinates for an equidistant grid,
1474  ! but we can temporarily also use them for getting indices
1475  ! in stretched grids
1476  xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
1477  xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;
1478 
1479  if(stagger_grid.and.phyboundblock(igrid).and.bcphys) then
1480  block=>psc(igrid)
1481  do idims=1,ndim
1482  ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
1483  ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
1484  {^ifthreed
1485  ! avoid using undetermined ghost cells at physical boundary edges
1486  if(idims == 1) then
1487  if(neighbor_type(-1,0,0,igrid)==neighbor_boundary .or. &
1488  neighbor_type(1,0,0,igrid)==neighbor_boundary) then
1489  if(neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixcomin2=ixcommin2
1490  if(neighbor_type(0,0,-1,igrid)==neighbor_boundary) ixcomin3=ixcommin3
1491  if(neighbor_type(0,1,0,igrid)==neighbor_boundary) ixcomax2=ixcommax2
1492  if(neighbor_type(0,0,1,igrid)==neighbor_boundary) ixcomax3=ixcommax3
1493  end if
1494  else if(idims == 2) then
1495  if(neighbor_type(0,-1,0,igrid)==neighbor_boundary .or. &
1496  neighbor_type(0,1,0,igrid)==neighbor_boundary) then
1497  if(neighbor_type(0,0,-1,igrid)==neighbor_boundary) ixcomin3=ixcommin3
1498  if(neighbor_type(0,0,1,igrid)==neighbor_boundary) ixcomax3=ixcommax3
1499  end if
1500  end if
1501  }
1502  do iside=1,2
1503  ii^d=kr(^d,idims)*(2*iside-3);
1504  if(neighbor_type(ii^d,igrid)/=neighbor_boundary) cycle
1505  if(( {(iside==1.and.idims==^d.and.ixcomin^d<ixcogmin^d+nghostcells)|.or.} ) &
1506  .or.( {(iside==2.and.idims==^d.and.ixcomax^d>ixcogmax^d-nghostcells)|.or. })) then
1507  {ixbmin^d=merge(ixcogmin^d,ixcomin^d,idims==^d);}
1508  {ixbmax^d=merge(ixcogmax^d,ixcomax^d,idims==^d);}
1509  call bc_phys(iside,idims,time,0.d0,psc(igrid),ixcog^l,ixb^l)
1510  end if
1511  end do
1512  end do
1513  end if
1514 
1515  if(prolongprimitive) then
1516  ! following line again assumes equidistant grid, but
1517  ! just computes indices, so also ok for stretched case
1518  ! reason for +1-1 and +1+1: the coarse representation has
1519  ! also nghostcells at each side. During
1520  ! prolongation, we need cells to left and right, hence -1/+1
1521  block=>psc(igrid)
1522  ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
1523  ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
1524  call phys_to_primitive(ixcog^l,ixco^l,psc(igrid)%w,psc(igrid)%x)
1525  end if
1526 
1527  if(ghost_copy) then
1528  call interpolation_copy(igrid,ixfi^l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
1529  else
1530  call interpolation_linear(igrid,ixfi^l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
1531  end if
1532 
1533  if(prolongprimitive) then
1534  block=>psc(igrid)
1535  call phys_to_conserved(ixcog^l,ixco^l,psc(igrid)%w,psc(igrid)%x)
1536  end if
1537 
1538  end subroutine bc_prolong
1539 
1540  subroutine bc_prolong_stg(igrid,i^D,iib^D,NeedProlong)
1541  use mod_amr_fct
1542  integer :: igrid,i^D,iib^D
1543  logical,dimension(-1:1^D&) :: NeedProlong
1544  logical :: fine_^Lin
1545  integer :: ixFi^L,ixCo^L
1546  double precision :: dxFi^D,dxCo^D,xFimin^D,xComin^D,invdxCo^D
1547  ! Check what is already at the desired level
1548  fine_^lin=.false.;
1549  {
1550  if(i^d>-1) fine_min^din=(.not.needprolong(i^dd-kr(^d,^dd)).and.neighbor_type(i^dd-kr(^d,^dd),igrid)/=1)
1551  if(i^d<1) fine_max^din=(.not.needprolong(i^dd+kr(^d,^dd)).and.neighbor_type(i^dd+kr(^d,^dd),igrid)/=1)
1552  \}
1553 
1554  ixfi^l=ixr_srl_^l(iib^d,i^d);
1555 
1556  dxfi^d=rnode(rpdx^d_,igrid);
1557  dxco^d=two*dxfi^d;
1558  invdxco^d=1.d0/dxco^d;
1559 
1560  xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
1561  xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;
1562 
1563  ! moved the physical boundary filling here, to only fill the
1564  ! part needed
1565 
1566  ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
1567  ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
1568 
1569  if(prolongprimitive) call phys_to_primitive(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)
1570 
1571  call prolong_2nd_stg(psc(igrid),psb(igrid),ixco^l,ixfi^l,dxco^d,xcomin^d,dxfi^d,xfimin^d,.true.,fine_^lin)
1572 
1573  if(prolongprimitive) call phys_to_conserved(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)
1574 
1575  ! The current region has already been refined, so it does not need to be prolonged again
1576  needprolong(i^d)=.false.
1577 
1578  end subroutine bc_prolong_stg
1579 
1580  subroutine interpolation_linear(igrid,ixFi^L,dxFi^D,xFimin^D, &
1581  dxCo^D,invdxCo^D,xComin^D)
1582  use mod_physics, only: phys_to_conserved
1583  integer, intent(in) :: igrid, ixFi^L
1584  double precision, intent(in) :: dxFi^D, xFimin^D,dxCo^D, invdxCo^D, xComin^D
1585 
1586  integer :: ixCo^D, jxCo^D, hxCo^D, ixFi^D, ix^D, iw, idims, nwmin,nwmax
1587  double precision :: xCo^D, xFi^D, eta^D
1588  double precision :: slopeL, slopeR, slopeC, signC, signR
1589  double precision :: slope(1:nw,ndim)
1590  !!double precision :: local_invdxCo^D
1591  double precision :: signedfactorhalf^D
1592  !integer :: ixshift^D, icase
1593 
1594  !icase=mod(nghostcells,2)
1595 
1596  if(prolongprimitive) then
1597  nwmin=1
1598  nwmax=nw
1599  else
1600  nwmin=nwhead
1601  nwmax=nwtail
1602  end if
1603 
1604  {do ixfi^db = ixfi^lim^db
1605  ! cell-centered coordinates of fine grid point
1606  ! here we temporarily use an equidistant grid
1607  xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
1608 
1609  ! indices of coarse cell which contains the fine cell
1610  ! since we computed lower left corner earlier
1611  ! in equidistant fashion: also ok for stretched case
1612  ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1
1613 
1614  ! cell-centered coordinates of coarse grid point
1615  ! here we temporarily use an equidistant grid
1616  xco^db=xcomin^db+(dble(ixco^db)-half)*dxco^db \}
1617 
1618  !if(.not.slab) then
1619  ! ^D&local_invdxCo^D=1.d0/psc(igrid)%dx({ixCo^DD},^D)\
1620  !endif
1621 
1622  if(slab_uniform) then
1623  ! actual cell-centered coordinates of fine grid point
1624  !!^D&xFi^D=block%x({ixFi^DD},^D)\
1625  ! actual cell-centered coordinates of coarse grid point
1626  !!^D&xCo^D=psc(igrid)%x({ixCo^DD},^D)\
1627  ! normalized distance between fine/coarse cell center
1628  ! in coarse cell: ranges from -0.5 to 0.5 in each direction
1629  ! (origin is coarse cell center)
1630  ! this is essentially +1/4 or -1/4 on cartesian mesh
1631  eta^d=(xfi^d-xco^d)*invdxco^d;
1632  else
1633  !select case(icase)
1634  ! case(0)
1635  !{! here we assume an even number of ghostcells!!!
1636  !ixshift^D=2*(mod(ixFi^D,2)-1)+1
1637  !if(ixshift^D>0.0d0)then
1638  ! ! oneven fine grid points
1639  ! eta^D=-0.5d0*(one-block%dvolume(ixFi^DD) &
1640  ! /sum(block%dvolume(ixFi^D:ixFi^D+1^D%ixFi^DD)))
1641  !else
1642  ! ! even fine grid points
1643  ! eta^D=+0.5d0*(one-block%dvolume(ixFi^DD) &
1644  ! /sum(block%dvolume(ixFi^D-1:ixFi^D^D%ixFi^DD)))
1645  !endif\}
1646  ! case(1)
1647  !{! here we assume an odd number of ghostcells!!!
1648  !ixshift^D=2*(mod(ixFi^D,2)-1)+1
1649  !if(ixshift^D>0.0d0)then
1650  ! ! oneven fine grid points
1651  ! eta^D=+0.5d0*(one-block%dvolume(ixFi^DD) &
1652  ! /sum(block%dvolume(ixFi^D-1:ixFi^D^D%ixFi^DD)))
1653  !else
1654  ! ! even fine grid points
1655  ! eta^D=-0.5d0*(one-block%dvolume(ixFi^DD) &
1656  ! /sum(block%dvolume(ixFi^D:ixFi^D+1^D%ixFi^DD)))
1657  !endif\}
1658  ! case default
1659  ! call mpistop("no such case")
1660  !end select
1661  ! the different cases for even/uneven number of ghost cells
1662  ! are automatically handled using the relative index to ixMlo
1663  ! as well as the pseudo-coordinates xFi and xCo
1664  ! these latter differ from actual cell centers when stretching is used
1665  ix^d=2*int((ixfi^d+ixmlo^d)/2)-ixmlo^d;
1666  {if(xfi^d>xco^d) then
1667  signedfactorhalf^d=0.5d0
1668  else
1669  signedfactorhalf^d=-0.5d0
1670  end if
1671  eta^d=signedfactorhalf^d*(one-psb(igrid)%dvolume(ixfi^dd) &
1672  /sum(psb(igrid)%dvolume(ix^d:ix^d+1^d%ixFi^dd))) \}
1673  !{eta^D=(xFi^D-xCo^D)*invdxCo^D &
1674  ! *two*(one-block%dvolume(ixFi^DD) &
1675  ! /sum(block%dvolume(ix^D:ix^D+1^D%ixFi^DD))) \}
1676  end if
1677 
1678  do idims=1,ndim
1679  hxco^d=ixco^d-kr(^d,idims)\
1680  jxco^d=ixco^d+kr(^d,idims)\
1681 
1682  do iw=nwmin,nwmax
1683  slopel=psc(igrid)%w(ixco^d,iw)-psc(igrid)%w(hxco^d,iw)
1684  sloper=psc(igrid)%w(jxco^d,iw)-psc(igrid)%w(ixco^d,iw)
1685  slopec=half*(sloper+slopel)
1686 
1687  ! get limited slope
1688  signr=sign(one,sloper)
1689  signc=sign(one,slopec)
1690  !select case(prolong_limiter)
1691  !case(1)
1692  ! ! unlimit
1693  ! slope(iw,idims)=slopeC
1694  !case(2)
1695  ! ! minmod
1696  ! slope(iw,idims)=signR*max(zero,min(dabs(slopeR), &
1697  ! signR*slopeL))
1698  !case(3)
1699  ! ! woodward
1700  ! slope(iw,idims)=two*signR*max(zero,min(dabs(slopeR), &
1701  ! signR*slopeL,signR*half*slopeC))
1702  !case(4)
1703  ! ! koren
1704  ! slope(iw,idims)=signR*max(zero,min(two*signR*slopeL, &
1705  ! (dabs(slopeR)+two*slopeL*signR)*third,two*dabs(slopeR)))
1706  !case default
1707  slope(iw,idims)=signc*max(zero,min(dabs(slopec), &
1708  signc*slopel,signc*sloper))
1709  !end select
1710  end do
1711  end do
1712 
1713  ! Interpolate from coarse cell using limited slopes
1714  psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)+&
1715  {(slope(nwmin:nwmax,^d)*eta^d)+}
1716 
1717  {end do\}
1718 
1719  if(prolongprimitive) then
1720  block=>psb(igrid)
1721  call phys_to_conserved(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)
1722  end if
1723 
1724  end subroutine interpolation_linear
1725 
1726  subroutine interpolation_copy(igrid, ixFi^L,dxFi^D,xFimin^D, &
1727  dxCo^D,invdxCo^D,xComin^D)
1728  use mod_physics, only: phys_to_conserved
1729  integer, intent(in) :: igrid, ixFi^L
1730  double precision, intent(in) :: dxFi^D, xFimin^D,dxCo^D, invdxCo^D, xComin^D
1731 
1732  integer :: ixCo^D, ixFi^D, nwmin,nwmax
1733  double precision :: xFi^D
1734 
1735  if(prolongprimitive) then
1736  nwmin=1
1737  nwmax=nw
1738  else
1739  nwmin=nwhead
1740  nwmax=nwtail
1741  end if
1742 
1743  {do ixfi^db = ixfi^lim^db
1744  ! cell-centered coordinates of fine grid point
1745  xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
1746 
1747  ! indices of coarse cell which contains the fine cell
1748  ! note: this also works for stretched grids
1749  ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1\}
1750 
1751  ! Copy from coarse cell
1752  psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)
1753 
1754  {end do\}
1755 
1756  if(prolongprimitive) call phys_to_conserved(ixg^ll,ixfi^l,psb(igrid)%w,psb(igrid)%x)
1757 
1758  end subroutine interpolation_copy
1759 
1760  subroutine pole_copy(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,ipole)
1761 
1762  integer, intent(in) :: ixIR^L,ixR^L,ixIS^L,ixS^L,ipole
1763  double precision :: wrecv(ixIR^S,1:nw), wsend(ixIS^S,1:nw)
1764 
1765  integer :: iw, iside, iB
1766 
1767  select case (ipole)
1768  {case (^d)
1769  iside=int((i^d+3)/2)
1770  ib=2*(^d-1)+iside
1771  do iw=nwhead,nwtail
1772  select case (typeboundary(iw,ib))
1773  case (bc_symm)
1774  wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1775  case (bc_asymm)
1776  wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1777  case default
1778  call mpistop("Pole boundary condition should be symm or asymm")
1779  end select
1780  end do \}
1781  end select
1782 
1783  end subroutine pole_copy
1784 
1785  subroutine pole_copy_stg(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,idirs,ipole)
1786 
1787  integer, intent(in) :: ixIR^L,ixR^L,ixIS^L,ixS^L,idirs,ipole
1788 
1789  double precision :: wrecv(ixIR^S,1:nws), wsend(ixIS^S,1:nws)
1790  integer :: iB, iside
1791 
1792  select case (ipole)
1793  {case (^d)
1794  iside=int((i^d+3)/2)
1795  ib=2*(^d-1)+iside
1796  select case (typeboundary(iw_mag(idirs),ib))
1797  case (bc_symm)
1798  wrecv(ixr^s,idirs) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
1799  case (bc_asymm)
1800  wrecv(ixr^s,idirs) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
1801  case default
1802  call mpistop("Pole boundary condition should be symm or asymm")
1803  end select
1804  \}
1805  end select
1806 
1807  end subroutine pole_copy_stg
1808 
1809  subroutine pole_buffer(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L)
1810 
1811  integer, intent(in) :: ixIR^L,ixR^L,ixIS^L,ixS^L
1812  double precision :: wrecv(ixIR^S,nwhead:nwtail), wsend(ixIS^S,1:nw)
1813 
1814  integer :: iw, iside, iB
1815 
1816  select case (ipole)
1817  {case (^d)
1818  iside=int((i^d+3)/2)
1819  ib=2*(^d-1)+iside
1820  do iw=nwhead,nwtail
1821  select case (typeboundary(iw,ib))
1822  case (bc_symm)
1823  wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1824  case (bc_asymm)
1825  wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1826  case default
1827  call mpistop("Pole boundary condition should be symm or asymm")
1828  end select
1829  end do \}
1830  end select
1831 
1832  end subroutine pole_buffer
1833 
1834  end subroutine getbc
1835 
1836  subroutine identifyphysbound(s,iib^D)
1838 
1839  type(state) :: s
1840  integer, intent(out) :: iib^D
1841 
1842  {
1843  if(s%is_physical_boundary(2*^d) .and. &
1844  s%is_physical_boundary(2*^d-1)) then
1845  iib^d=2
1846  else if(s%is_physical_boundary(2*^d-1)) then
1847  iib^d=-1
1848  else if(s%is_physical_boundary(2*^d)) then
1849  iib^d=1
1850  else
1851  iib^d=0
1852  end if
1853  \}
1854 
1855  end subroutine identifyphysbound
1856 
1857 end module mod_ghostcells_update
subroutine bc_phys(iside, idims, time, qdt, s, ixGL, ixBL)
fill ghost cells at a physical boundary
subroutine getintbc(time, ixGL)
fill inner boundary values
subroutine coarsen_grid(sFi, ixFiGL, ixFiL, sCo, ixCoGL, ixCoL)
coarsen one grid to its coarser representative
Definition: coarsen.t:3
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
subroutine bc_recv_restrict
Receive from fine neighbor.
subroutine gc_prolong(igrid)
subroutine interpolation_linear(igrid, ixFiL, dxFiD, xFiminD, dxCoD, invdxCoD, xCominD)
subroutine bc_prolong(igrid, iD, iibD)
do prolongation for fine blocks after receipt data from coarse neighbors
subroutine indices_for_syncing(idir, iD, ixRL, ixSL, ixRsyncL, ixSsyncL)
subroutine bc_fill_srl_stg
fill siblings ghost cells with received data
subroutine pole_copy(wrecv, ixIRL, ixRL, wsend, ixISL, ixSL, ipole)
subroutine bc_recv_prolong
Receive from coarse neighbor.
subroutine bc_fill_prolong(igrid, iD, iibD)
Send to finer neighbor.
subroutine bc_recv_srl
Receive from sibling at same refinement level.
subroutine bc_fill_restrict(igrid, iD, iibD)
fill coarser neighbor's ghost cells
subroutine fill_coarse_boundary(igrid, iD)
subroutine interpolation_copy(igrid, ixFiL, dxFiD, xFiminD, dxCoD, invdxCoD, xCominD)
subroutine pole_copy_stg(wrecv, ixIRL, ixRL, wsend, ixISL, ixSL, idirs, ipole)
subroutine bc_send_restrict
Send to coarser neighbor.
subroutine pole_buffer(wrecv, ixIRL, ixRL, wsend, ixISL, ixSL)
subroutine fill_boundary_before_gc(igrid)
Physical boundary conditions.
subroutine bc_fill_srl(igrid, iD, iibD)
subroutine bc_fill_restrict_stg
fill restricted ghost cells after receipt
subroutine bc_send_srl
Send to sibling at same refinement level.
subroutine fill_boundary_after_gc(igrid)
Physical boundary conditions.
logical function skip_direction(dir)
subroutine bc_send_prolong
Send to finer neighbor.
subroutine bc_fill_prolong_stg
fill coarser representative with data from coarser neighbors
subroutine bc_prolong_stg(igrid, iD, iibD, NeedProlong)
subroutine, public prolong_2nd_stg(sCo, sFi, ixCoLin, ixFiLin, dxCoD, xCominD, dxFiD, xFiminD, ghost, fine_Lin)
This subroutine performs a 2nd order prolongation for a staggered field F, preserving the divergence ...
Definition: mod_amr_fct.t:41
update ghost cells of all blocks including physical boundaries
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &) sizes_r_send_total
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p2
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p2
integer, dimension(^nd, 0:3) l
integer, dimension(:), allocatable sendrequest_c_p
integer, dimension(^nd,-1:1) ixs_r_stg_
integer, dimension(^nd,-1:1) ixs_srl_stg_
subroutine get_bc_comm_type(comm_type, ixL, ixGL, nwstart, nwbc)
integer, dimension(^nd, 0:3^d &) sizes_p_send_stg
integer, dimension(-1:1^d &) sizes_srl_send_total
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p1
subroutine identifyphysbound(s, iibD)
integer, dimension(^nd, 0:3) ixr_r_stg_
integer, dimension(:), allocatable sendrequest_r
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p1
integer, dimension(:), allocatable recvrequest_r
double precision, dimension(:), allocatable sendbuffer_p
integer, dimension(:,:), allocatable sendstatus_c_p
integer, dimension(^nd,-1:1^d &) sizes_r_send_stg
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_f
integer, dimension(-1:2,-1:1) ixr_srl_
integer, dimension(:^d &,:^d &), pointer type_recv_srl
integer, dimension(:,:), allocatable sendstatus_c_sr
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p2
double precision, dimension(:), allocatable recvbuffer_r
integer, dimension(:), allocatable recvrequest_c_p
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(:), allocatable recvrequest_c_sr
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p1
integer, dimension(:), allocatable sendrequest_c_sr
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p2
integer, dimension(-1:2,-1:1) ixs_srl_
double precision, dimension(:), allocatable recvbuffer_srl
integer, dimension(:^d &,:^d &), pointer type_send_srl
integer, dimension(:,:), allocatable recvstatus_p
integer, dimension(^nd, 0:3) ixs_p_stg_
integer, dimension(^nd, 0:3^d &) sizes_p_recv_stg
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p2
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(^nd,-1:1^d &) sizes_srl_recv_stg
double precision, dimension(:), allocatable recvbuffer_p
integer, dimension(^nd,-1:1) ixr_srl_stg_
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p1
double precision, dimension(:), allocatable sendbuffer_r
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_f
integer, dimension(-1:1, 0:3) ixr_p_
integer, dimension(-1:1, 0:3) ixr_r_
integer, dimension(0:3^d &) sizes_p_recv_total
integer, dimension(:,:), allocatable recvstatus_c_sr
integer, dimension(0:3^d &) sizes_r_recv_total
integer, parameter npwbuf
integer, dimension(-1:1, 0:3) ixs_p_
integer, dimension(:), allocatable sendrequest_p
integer, dimension(^nd, 0:3^d &) sizes_r_recv_stg
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
double precision, dimension(:), allocatable sendbuffer_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p2
integer, dimension(:,:), allocatable sendstatus_srl
integer, dimension(:,:), allocatable recvstatus_srl
integer, dimension(:), allocatable sendrequest_srl
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_f
integer, dimension(:,:), allocatable recvstatus_r
integer, dimension(-1:1,-1:1) ixs_r_
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p1
integer, dimension(0:3^d &) sizes_p_send_total
integer, dimension(:^d &,:^d &), pointer type_recv_p
integer, dimension(^nd, 0:3) ixr_p_stg_
integer, dimension(:^d &,:^d &), pointer type_send_p
integer, dimension(:,:), allocatable sendstatus_r
integer, dimension(:,:), allocatable sendstatus_p
integer, dimension(-1:1^d &) sizes_srl_recv_total
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
integer, dimension(:), allocatable recvrequest_p
integer, dimension(:), allocatable recvrequest_srl
integer, dimension(^nd,-1:1^d &) sizes_srl_send_stg
integer, dimension(:,:), allocatable recvstatus_c_p
This module contains definitions of global parameters and variables and some generic functions/subrou...
logical internalboundary
if there is an internal boundary
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
logical, dimension(:), allocatable phyboundblock
True if a block has any physical boundary.
integer, dimension(:), allocatable, parameter d
logical ghost_copy
whether copy values instead of interpolation in ghost cells of finer blocks
integer ierrmpi
A global MPI error return code.
integer nghostcells
Number of ghost cells surrounding a grid.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:59
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:27
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:58
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:19