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