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