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