121 integer :: nghostcellsCo, interpolation_order
122 integer :: nx^D, nxCo^D, ixG^L, i^D, idir
125 ixm^l=ixg^l^lsubnghostcells;
129 ixcogsmax^d=ixcogmax^d;
133 nx^d=ixmmax^d-ixmmin^d+1;
137 interpolation_order=1
139 interpolation_order=2
143 if (nghostcellsco+interpolation_order-1>
nghostcells)
then
144 call mpistop(
"interpolation order for prolongation in getbc too high")
150 ixs_srl_min^d(-1)=ixmmin^d
151 ixs_srl_min^d( 0)=ixmmin^d
154 ixs_srl_max^d( 0)=ixmmax^d
155 ixs_srl_max^d( 1)=ixmmax^d
158 ixr_srl_min^d( 0)=ixmmin^d
159 ixr_srl_min^d( 1)=ixmmax^d+1
161 ixr_srl_max^d( 0)=ixmmax^d
162 ixr_srl_max^d( 1)=ixgmax^d
164 ixs_r_min^d(-1)=ixcommin^d
165 ixs_r_min^d( 0)=ixcommin^d
168 ixs_r_max^d( 0)=ixcommax^d
169 ixs_r_max^d( 1)=ixcommax^d
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
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
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)
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)
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)
201 allocate(pole_buf%ws(ixgs^t,nws))
204 { ixs_srl_stg_min^d(idir,-1)=ixmmin^d-
kr(idir,^d)
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
211 ixr_srl_stg_min^d(idir,-1)=1-
kr(idir,^d)
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
218 ixs_r_stg_min^d(idir,-1)=ixcommin^d-
kr(idir,^d)
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
225 ixr_r_stg_min^d(idir,0)=1-
kr(idir,^d)
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
237 ixs_p_stg_min^d(idir,0)=ixmmin^d-1
238 ixs_p_stg_max^d(idir,0)=ixmmin^d-1+nghostcellsco
239 ixs_p_stg_min^d(idir,1)=ixmmin^d-1
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
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
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
253 ixr_p_stg_max^d(idir,3)=ixcommax^d+nghostcellsco
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
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)
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)|*}
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)|*}
350 subroutine getbc(time,qdt,psb,nwstart,nwbc)
357 double precision,
intent(in) :: time, qdt
358 type(state),
target :: psb(max_blocks)
359 integer,
intent(in) :: nwstart
360 integer,
intent(in) :: nwbc
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
369 integer :: ibuf_start, ibuf_next
371 integer,
dimension(1) :: shapes
374 time_bcin=mpi_wtime()
377 nwtail=nwstart+nwbc-1
386 do iigrid=1,igridstail; igrid=igrids(iigrid);
387 if(any(neighbor_type(:^d&,igrid)==neighbor_coarse))
then
416 do iigrid=1,igridstail; igrid=igrids(iigrid);
419 select case (neighbor_type(i^d,igrid))
420 case (neighbor_sibling)
423 call bc_recv_restrict
429 do iigrid=1,igridstail; igrid=igrids(iigrid);
432 select case (neighbor_type(i^d,igrid))
433 case (neighbor_sibling)
435 case (neighbor_coarse)
436 call bc_send_restrict
443 do iigrid=1,igridstail; igrid=igrids(iigrid);
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)
459 if(stagger_grid)
then
467 do iigrid=1,igridstail; igrid=igrids(iigrid);
470 select case (neighbor_type(i^d,igrid))
471 case (neighbor_sibling)
474 call bc_fill_restrict_stg
481 if (isend_buf(ipwbuf)/=0)
deallocate(pwbuf(ipwbuf)%w)
490 do iigrid=1,igridstail; igrid=igrids(iigrid);
493 if (neighbor_type(i^d,igrid)==neighbor_coarse)
call bc_recv_prolong
497 do iigrid=1,igridstail; igrid=igrids(iigrid);
500 if (neighbor_type(i^d,igrid)==neighbor_fine)
call bc_send_prolong
506 do iigrid=1,igridstail; igrid=igrids(iigrid);
509 if (neighbor_type(i^d,igrid)==neighbor_fine)
call bc_fill_prolong(igrid,i^d)
517 if(stagger_grid)
then
523 do iigrid=1,igridstail; igrid=igrids(iigrid);
526 if(neighbor_type(i^d,igrid)==neighbor_coarse)
call bc_fill_prolong_stg
532 do iigrid=1,igridstail; igrid=igrids(iigrid);
533 call gc_prolong(igrid)
538 if (isend_buf(ipwbuf)/=0)
deallocate(pwbuf(ipwbuf)%w)
544 do iigrid=1,igridstail; igrid=igrids(iigrid);
545 if(.not.phyboundblock(igrid)) cycle
546 call fill_boundary_after_gc(igrid)
552 if(
bcphys.and.
associated(phys_boundary_adjust))
then
554 do iigrid=1,igridstail; igrid=igrids(iigrid);
555 if(.not.phyboundblock(igrid)) cycle
556 call phys_boundary_adjust(igrid,psb)
561 time_bc=time_bc+(mpi_wtime()-time_bcin)
566 integer,
intent(in) :: dir(^nd)
568 if (all(dir == 0))
then
576 subroutine fill_boundary_after_gc(igrid)
578 integer,
intent(in) :: igrid
580 integer :: idims,iside,i^d,k^
l,ixb^
l
583 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
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)}
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;
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
605 if (neighbor_type(i^d,igrid) /= neighbor_boundary) cycle
607 call bc_phys(iside,idims,time,qdt,psb(igrid),ixg^ll,ixb^
l,nwhead,nwtail)
611 end subroutine fill_boundary_after_gc
614 subroutine bc_recv_srl
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)+}
622 if(stagger_grid)
then
630 end subroutine bc_recv_srl
633 subroutine bc_recv_restrict
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), &
643 if(stagger_grid)
then
646 mpi_double_precision,ipe_neighbor,itag, &
653 end subroutine bc_recv_restrict
656 subroutine bc_send_srl
658 ipe_neighbor=neighbor(2,i^d,igrid)
660 if(ipe_neighbor/=mype)
then
661 ineighbor=neighbor(1,i^d,igrid)
662 ipole=neighbor_pole(i^d,igrid)
666 itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
669 if(stagger_grid)
then
676 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
688 n_i^d=i^d^d%n_i^dd=-i^dd;\}
690 if (isend_buf(ipwbuf)/=0)
then
693 deallocate(pwbuf(ipwbuf)%w)
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)
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, &
703 ipwbuf=1+modulo(ipwbuf,
npwbuf)
704 if(stagger_grid)
then
711 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
722 end subroutine bc_send_srl
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
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)
736 psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
737 psb(igrid)%w(ixs^s,nwhead:nwtail)
738 if(stagger_grid)
then
742 psb(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
749 n_i^d=i^d^d%n_i^dd=-i^dd;\}
752 call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^
l,psb(igrid)%w,ixg^ll,ixs^
l,ipole)
753 if(stagger_grid)
then
757 call pole_copy_stg(psb(ineighbor)%ws,ixgs^ll,ixr^
l,psb(igrid)%ws,ixgs^ll,ixs^
l,idir,ipole)
763 end subroutine bc_fill_srl
766 subroutine bc_send_restrict
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)
777 itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
780 if(stagger_grid)
then
787 reshape(psc(igrid)%ws(ixs^s,idir),shapes)
792 mpi_double_precision,ipe_neighbor,itag, &
800 n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
802 if(isend_buf(ipwbuf)/=0)
then
805 deallocate(pwbuf(ipwbuf)%w)
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)
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, &
815 ipwbuf=1+modulo(ipwbuf,
npwbuf)
816 if(stagger_grid)
then
823 reshape(psc(igrid)%ws(ixs^s,idir),shapes)
828 mpi_double_precision,ipe_neighbor,itag, &
835 end subroutine bc_send_restrict
838 subroutine bc_fill_restrict(igrid,i^D)
839 integer,
intent(in) :: igrid,i^d
841 integer :: ic^d,n_inc^d,ixs^
l,ixr^
l,ipe_neighbor,ineighbor,ipole,idir
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)
853 psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
854 psc(igrid)%w(ixs^s,nwhead:nwtail)
855 if(stagger_grid)
then
859 psb(ineighbor)%ws(ixr^s,idir)=psc(igrid)%ws(ixs^s,idir)
866 n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
869 call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^
l,psc(igrid)%w,
ixcog^
l,ixs^
l,ipole)
870 if(stagger_grid)
then
875 call pole_copy_stg(psb(ineighbor)%ws,ixgs^ll,ixr^
l,psc(igrid)%ws,
ixcogs^
l,ixs^
l,idir,ipole)
881 end subroutine bc_fill_restrict
884 subroutine bc_fill_srl_stg
885 integer :: ixs^
l,ixr^
l,n_i^d,ixssync^
l,ixrsync^
l,idir
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)
902 shape=shape(psb(igrid)%ws(ixs^s,idir)))
908 n_i^d=i^d^d%n_i^dd=-i^dd;\}
916 shape=shape(psb(igrid)%ws(ixs^s,idir)))
918 call pole_copy_stg(psb(igrid)%ws,ixgs^ll,ixr^
l,pole_buf%ws,ixgs^ll,ixs^
l,idir,ipole)
923 end subroutine bc_fill_srl_stg
926 subroutine bc_fill_restrict_stg
928 ipole=neighbor_pole(i^d,igrid)
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)
942 shape=shape(psb(igrid)%ws(ixr^s,idir)))
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)
955 n_i^d=i^d^d%n_i^dd=-i^dd;\}
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)
973 end subroutine bc_fill_restrict_stg
976 subroutine bc_recv_prolong
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
981 ipe_neighbor=neighbor(2,i^d,igrid)
982 if (ipe_neighbor/=mype)
then
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), &
988 if(stagger_grid)
then
991 mpi_double_precision,ipe_neighbor,itag,&
997 end subroutine bc_recv_prolong
1000 subroutine bc_send_prolong
1003 ipole=neighbor_pole(i^d,igrid)
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)
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), &
1017 if(stagger_grid)
then
1024 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
1025 ibuf_start=ibuf_next
1029 mpi_double_precision,ipe_neighbor,itag, &
1037 n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
1039 if(isend_buf(ipwbuf)/=0)
then
1042 deallocate(pwbuf(ipwbuf)%w)
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)
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, &
1052 ipwbuf=1+modulo(ipwbuf,
npwbuf)
1053 if(stagger_grid)
then
1060 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
1061 ibuf_start=ibuf_next
1065 mpi_double_precision,ipe_neighbor,itag, &
1073 end subroutine bc_send_prolong
1076 subroutine bc_fill_prolong(igrid,i^D)
1077 integer,
intent(in) :: igrid,i^d
1079 integer :: ipe_neighbor,ineighbor,ixs^
l,ixr^
l,ic^d,inc^d,ipole,idir
1081 ipole=neighbor_pole(i^d,igrid)
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
1089 ineighbor=neighbor_child(1,inc^d,igrid)
1093 psc(ineighbor)%w(ixr^s,nwhead:nwtail) &
1094 =psb(igrid)%w(ixs^s,nwhead:nwtail)
1095 if(stagger_grid)
then
1099 psc(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
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
1110 ineighbor=neighbor_child(1,inc^d,igrid)
1113 n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
1116 call pole_copy(psc(ineighbor)%w,
ixcog^
l,ixr^
l,psb(igrid)%w,ixg^ll,ixs^
l,ipole)
1117 if(stagger_grid)
then
1121 call pole_copy_stg(psc(ineighbor)%ws,
ixcogs^
l,ixr^
l,psb(igrid)%ws,ixgs^ll,ixs^
l,idir,ipole)
1127 end subroutine bc_fill_prolong
1129 subroutine gc_prolong(igrid)
1130 integer,
intent(in) :: igrid
1132 integer :: i^d,idims,iside
1133 logical,
dimension(-1:1^D&) :: needprolong
1138 if (neighbor_type(i^d,igrid)==neighbor_coarse)
then
1139 call bc_prolong(igrid,i^d)
1140 needprolong(i^d)=.true.
1143 if(stagger_grid)
then
1154 if (needprolong(i^dd))
call bc_prolong_stg(igrid,i^dd,needprolong)
1165 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1171 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1177 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1183 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1186 end subroutine gc_prolong
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
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)
1205 shape=shape(psc(igrid)%ws(ixr^s,idir)))
1212 n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
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)
1227 end subroutine bc_fill_prolong_stg
1230 subroutine bc_prolong(igrid,i^D)
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
1238 dxfi^d=rnode(rpdx^d_,igrid);
1240 invdxco^d=1.d0/dxco^d;
1246 xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
1247 xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;
1249 if(phyboundblock(igrid).and.
bcphys)
then
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;
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
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
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)
1285 if(prolongprimitive)
then
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;
1298 call interpolation_copy(igrid,ixfi^
l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
1300 call interpolation_linear(igrid,ixfi^
l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
1303 if(prolongprimitive)
then
1308 end subroutine bc_prolong
1310 subroutine bc_prolong_stg(igrid,i^D,NeedProlong)
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
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)
1326 dxfi^d=rnode(rpdx^d_,igrid);
1328 invdxco^d=1.d0/dxco^d;
1330 xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
1331 xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;
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;
1339 if(prolongprimitive)
call phys_to_primitive(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1341 call prolong_2nd_stg(psc(igrid),psb(igrid),ixco^
l,ixfi^
l,dxco^d,xcomin^d,dxfi^d,xfimin^d,.true.,fine_^lin)
1343 if(prolongprimitive)
call phys_to_conserved(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1346 needprolong(i^d)=.false.
1348 end subroutine bc_prolong_stg
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
1356 double precision :: xco^d, xfi^d, eta^d
1357 double precision :: slopel, sloper, slopec, signc, signr
1358 double precision :: slope(1:nw,ndim)
1360 double precision :: signedfactorhalf^d
1361 integer :: ixco^d, jxco^d, hxco^d, ixfi^d, ix^d, iw, idims, nwmin,nwmax
1366 if(prolongprimitive)
then
1374 {
do ixfi^db = ixfi^lim^db
1377 xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
1382 ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1
1386 xco^db=xcomin^db+(dble(ixco^db)-half)*dxco^db \}
1392 if(slab_uniform)
then
1401 eta^d=(xfi^d-xco^d)*invdxco^d;
1435 ix^d=2*int((ixfi^d+ixmlo^d)/2)-ixmlo^d;
1436 {
if(xfi^d>xco^d)
then
1437 signedfactorhalf^d=0.5d0
1439 signedfactorhalf^d=-0.5d0
1441 eta^d=signedfactorhalf^d*(one-psb(igrid)%dvolume(ixfi^dd) &
1442 /sum(psb(igrid)%dvolume(ix^d:ix^d+1^d%ixFi^dd))) \}
1449 hxco^d=ixco^d-kr(^d,idims)\
1450 jxco^d=ixco^d+kr(^d,idims)\
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)
1458 signr=sign(one,sloper)
1459 signc=sign(one,slopec)
1477 slope(iw,idims)=signc*max(zero,min(dabs(slopec), &
1478 signc*slopel,signc*sloper))
1484 psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)+&
1485 {(slope(nwmin:nwmax,^d)*eta^d)+}
1489 if(prolongprimitive)
then
1491 call phys_to_conserved(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1494 end subroutine interpolation_linear
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
1502 double precision :: xfi^d
1503 integer :: ixco^d, ixfi^d, nwmin,nwmax
1505 if(prolongprimitive)
then
1513 {
do ixfi^db = ixfi^lim^db
1515 xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
1519 ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1\}
1522 psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)
1526 if(prolongprimitive)
call phys_to_conserved(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1528 end subroutine interpolation_copy
1530 subroutine pole_copy(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,ipole)
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)
1535 integer :: iw, iside, ib
1539 iside=int((i^d+3)/2)
1542 select case (typeboundary(iw,ib))
1544 wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1546 wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1548 call mpistop(
"Pole boundary condition should be symm or asymm")
1553 end subroutine pole_copy
1555 subroutine pole_copy_stg(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,idirs,ipole)
1557 integer,
intent(in) :: ixir^
l,ixr^
l,ixis^
l,ixs^
l,idirs,ipole
1559 double precision :: wrecv(ixir^s,1:nws), wsend(ixis^s,1:nws)
1560 integer :: ib, iside
1564 iside=int((i^d+3)/2)
1566 select case (typeboundary(iw_mag(idirs),ib))
1568 wrecv(ixr^s,idirs) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
1570 wrecv(ixr^s,idirs) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
1572 call mpistop(
"Pole boundary condition should be symm or asymm")
1577 end subroutine pole_copy_stg
1579 subroutine pole_buffer(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L)
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)
1584 integer :: iw, iside, ib
1588 iside=int((i^d+3)/2)
1591 select case (typeboundary(iw,ib))
1593 wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1595 wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1597 call mpistop(
"Pole boundary condition should be symm or asymm")
1602 end subroutine pole_buffer