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 :: nwhead, nwtail
364 integer :: iigrid, igrid, isizes, i^D
365 integer :: isend_buf(npwbuf), ipwbuf, nghostcellsco
367 integer :: ibuf_start, ibuf_next
369 integer,
dimension(1) :: shapes
372 time_bcin=mpi_wtime()
375 nwtail=nwstart+nwbc-1
384 do iigrid=1,igridstail; igrid=igrids(iigrid);
385 if(any(neighbor_type(:^d&,igrid)==neighbor_coarse))
then
415 do iigrid=1,igridstail; igrid=igrids(iigrid);
418 select case (neighbor_type(i^d,igrid))
419 case (neighbor_sibling)
420 call bc_recv_srl(igrid,i^d)
422 call bc_recv_restrict(igrid,i^d)
428 do iigrid=1,igridstail; igrid=igrids(iigrid);
431 select case (neighbor_type(i^d,igrid))
432 case (neighbor_sibling)
433 call bc_send_srl(igrid,i^d)
434 case (neighbor_coarse)
435 call bc_send_restrict(igrid,i^d)
442 if(stagger_grid)
then
450 if (isend_buf(ipwbuf)/=0)
deallocate(pwbuf(ipwbuf)%w)
454 do iigrid=1,igridstail; igrid=igrids(iigrid);
457 select case (neighbor_type(i^d,igrid))
458 case(neighbor_sibling)
459 call bc_fill_srl(igrid,i^d)
460 case(neighbor_coarse)
461 call bc_fill_restrict(igrid,i^d)
467 if(stagger_grid)
then
471 do iigrid=1,igridstail; igrid=igrids(iigrid);
474 select case (neighbor_type(i^d,igrid))
475 case (neighbor_sibling)
476 call bc_fill_srl_stg(igrid,i^d)
478 call bc_fill_restrict_stg(igrid,i^d)
492 do iigrid=1,igridstail; igrid=igrids(iigrid);
495 if (neighbor_type(i^d,igrid)==neighbor_coarse)
call bc_recv_prolong(igrid,i^d)
499 do iigrid=1,igridstail; igrid=igrids(iigrid);
502 if (neighbor_type(i^d,igrid)==neighbor_fine)
call bc_send_prolong(igrid,i^d)
509 if(stagger_grid)
then
515 if (isend_buf(ipwbuf)/=0)
deallocate(pwbuf(ipwbuf)%w)
520 do iigrid=1,igridstail; igrid=igrids(iigrid);
523 if (neighbor_type(i^d,igrid)==neighbor_fine)
call bc_fill_prolong(igrid,i^d)
528 if(stagger_grid)
then
531 do iigrid=1,igridstail; igrid=igrids(iigrid);
534 if(neighbor_type(i^d,igrid)==neighbor_coarse)
call bc_fill_prolong_stg(igrid,i^d)
540 do iigrid=1,igridstail; igrid=igrids(iigrid);
541 call gc_prolong(igrid)
548 do iigrid=1,igridstail; igrid=igrids(iigrid);
549 if(.not.phyboundblock(igrid)) cycle
550 call fill_boundary_after_gc(igrid)
556 if(
bcphys.and.
associated(phys_boundary_adjust))
then
558 do iigrid=1,igridstail; igrid=igrids(iigrid);
559 if(.not.phyboundblock(igrid)) cycle
560 call phys_boundary_adjust(igrid,psb)
565 time_bc=time_bc+(mpi_wtime()-time_bcin)
570 integer,
intent(in) :: dir(^nd)
572 if (all(dir == 0))
then
580 subroutine fill_boundary_after_gc(igrid)
582 integer,
intent(in) :: igrid
584 integer :: idims,iside,i^d,k^
l,ixb^
l
587 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
594 kmin2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0,-1,igrid)==1)
595 kmax2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0, 1,igrid)==1)}
597 kmin2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0,-1,0,igrid)==1)
598 kmax2=merge(1, 0, idims .lt. 2 .and. neighbor_type(0, 1,0,igrid)==1)
599 kmin3=merge(1, 0, idims .lt. 3 .and. neighbor_type(0,0,-1,igrid)==1)
600 kmax3=merge(1, 0, idims .lt. 3 .and. neighbor_type(0,0, 1,igrid)==1)}
601 ixbmin^d=ixglo^d+kmin^d*nghostcells;
602 ixbmax^d=ixghi^d-kmax^d*nghostcells;
604 i^d=kr(^d,idims)*(2*iside-3);
605 if (aperiodb(idims))
then
606 if (neighbor_type(i^d,igrid) /= neighbor_boundary .and. &
607 .not. psb(igrid)%is_physical_boundary(2*idims-2+iside)) cycle
609 if (neighbor_type(i^d,igrid) /= neighbor_boundary) cycle
611 call bc_phys(iside,idims,time,qdt,psb(igrid),ixg^ll,ixb^
l)
615 end subroutine fill_boundary_after_gc
618 subroutine bc_recv_srl(igrid,i^D)
619 integer,
intent(in) :: igrid,i^d
621 integer :: ipe_neighbor
623 ipe_neighbor=neighbor(2,i^d,igrid)
624 if (ipe_neighbor/=mype)
then
626 itag=(3**^nd+4**^nd)*(igrid-1)+{(i^d+1)*3**(^d-1)+}
629 if(stagger_grid)
then
637 end subroutine bc_recv_srl
640 subroutine bc_recv_restrict(igrid,i^D)
641 integer,
intent(in) :: igrid,i^d
643 integer :: ic^d,inc^d,ipe_neighbor
645 {
do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
646 inc^db=2*i^db+ic^db\}
647 ipe_neighbor=neighbor_child(2,inc^d,igrid)
648 if (ipe_neighbor/=mype)
then
650 itag=(3**^nd+4**^nd)*(igrid-1)+3**^nd+{inc^d*4**(^d-1)+}
651 call mpi_irecv(psb(igrid)%w,1,
type_recv_r(inc^d), &
653 if(stagger_grid)
then
656 mpi_double_precision,ipe_neighbor,itag, &
663 end subroutine bc_recv_restrict
666 subroutine bc_send_srl(igrid,i^D)
667 integer,
intent(in) :: igrid,i^d
669 integer :: n_i^d,ipole,idir,ineighbor,ipe_neighbor,ixs^
l
671 ipe_neighbor=neighbor(2,i^d,igrid)
672 if(ipe_neighbor/=mype)
then
673 ineighbor=neighbor(1,i^d,igrid)
674 ipole=neighbor_pole(i^d,igrid)
678 itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
681 if(stagger_grid)
then
688 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
700 n_i^d=i^d^d%n_i^dd=-i^dd;\}
702 if (isend_buf(ipwbuf)/=0)
then
705 deallocate(pwbuf(ipwbuf)%w)
707 allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
708 call pole_buffer(pwbuf(ipwbuf)%w,ixs^
l,ixs^
l,psb(igrid)%w,ixg^ll,ixs^
l,ipole)
711 itag=(3**^nd+4**^nd)*(ineighbor-1)+{(n_i^d+1)*3**(^d-1)+}
712 isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
713 call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
715 ipwbuf=1+modulo(ipwbuf,
npwbuf)
716 if(stagger_grid)
then
723 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
734 end subroutine bc_send_srl
737 subroutine bc_send_restrict(igrid,i^D)
738 integer,
intent(in) :: igrid,i^d
740 integer :: ic^d,n_inc^d,ipole,idir,ineighbor,ipe_neighbor,ixs^
l
742 ipe_neighbor=neighbor(2,i^d,igrid)
743 if(ipe_neighbor/=mype)
then
744 ic^d=1+modulo(node(pig^d_,igrid)-1,2);
745 if({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.})
return
746 ineighbor=neighbor(1,i^d,igrid)
747 ipole=neighbor_pole(i^d,igrid)
751 itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
754 if(stagger_grid)
then
761 reshape(psc(igrid)%ws(ixs^s,idir),shapes)
766 mpi_double_precision,ipe_neighbor,itag, &
774 n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
776 if(isend_buf(ipwbuf)/=0)
then
779 deallocate(pwbuf(ipwbuf)%w)
781 allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
782 call pole_buffer(pwbuf(ipwbuf)%w,ixs^
l,ixs^
l,psc(igrid)%w,
ixcog^
l,ixs^
l,ipole)
785 itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
786 isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
787 call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
789 ipwbuf=1+modulo(ipwbuf,
npwbuf)
790 if(stagger_grid)
then
797 reshape(psc(igrid)%ws(ixs^s,idir),shapes)
802 mpi_double_precision,ipe_neighbor,itag, &
809 end subroutine bc_send_restrict
812 subroutine bc_fill_srl(igrid,i^D)
813 integer,
intent(in) :: igrid,i^d
815 integer :: ineighbor,ipe_neighbor,ipole,ixs^
l,ixr^
l,n_i^d,idir
817 ipe_neighbor=neighbor(2,i^d,igrid)
818 if(ipe_neighbor==mype)
then
819 ineighbor=neighbor(1,i^d,igrid)
820 ipole=neighbor_pole(i^d,igrid)
825 psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
826 psb(igrid)%w(ixs^s,nwhead:nwtail)
827 if(stagger_grid)
then
831 psb(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
838 n_i^d=i^d^d%n_i^dd=-i^dd;\}
841 call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^
l,psb(igrid)%w,ixg^ll,ixs^
l,ipole)
842 if(stagger_grid)
then
846 call pole_copy_stg(psb(ineighbor)%ws,ixgs^ll,ixr^
l,psb(igrid)%ws,ixgs^ll,ixs^
l,idir,ipole)
852 end subroutine bc_fill_srl
855 subroutine bc_fill_restrict(igrid,i^D)
856 integer,
intent(in) :: igrid,i^d
858 integer :: ic^d,n_inc^d,ixs^
l,ixr^
l,ipe_neighbor,ineighbor,ipole,idir
860 ipe_neighbor=neighbor(2,i^d,igrid)
861 if(ipe_neighbor==mype)
then
862 ic^d=1+modulo(node(pig^d_,igrid)-1,2);
863 if({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.})
return
864 ineighbor=neighbor(1,i^d,igrid)
865 ipole=neighbor_pole(i^d,igrid)
870 psb(ineighbor)%w(ixr^s,nwhead:nwtail)=&
871 psc(igrid)%w(ixs^s,nwhead:nwtail)
872 if(stagger_grid)
then
876 psb(ineighbor)%ws(ixr^s,idir)=psc(igrid)%ws(ixs^s,idir)
883 n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
886 call pole_copy(psb(ineighbor)%w,ixg^ll,ixr^
l,psc(igrid)%w,
ixcog^
l,ixs^
l,ipole)
887 if(stagger_grid)
then
892 call pole_copy_stg(psb(ineighbor)%ws,ixgs^ll,ixr^
l,psc(igrid)%ws,
ixcogs^
l,ixs^
l,idir,ipole)
898 end subroutine bc_fill_restrict
901 subroutine bc_fill_srl_stg(igrid,i^D)
902 integer,
intent(in) :: igrid,i^d
904 integer :: ixs^
l,ixr^
l,n_i^d,idir,ineighbor,ipe_neighbor,ipole
906 ipe_neighbor=neighbor(2,i^d,igrid)
907 if(ipe_neighbor/=mype)
then
908 ineighbor=neighbor(1,i^d,igrid)
909 ipole=neighbor_pole(i^d,igrid)
921 shape=shape(psb(igrid)%ws(ixs^s,idir)))
927 n_i^d=i^d^d%n_i^dd=-i^dd;\}
935 shape=shape(psb(igrid)%ws(ixs^s,idir)))
937 call pole_copy_stg(psb(igrid)%ws,ixgs^ll,ixr^
l,pole_buf%ws,ixgs^ll,ixs^
l,idir,ipole)
942 end subroutine bc_fill_srl_stg
945 subroutine bc_fill_restrict_stg(igrid,i^D)
946 integer,
intent(in) :: igrid,i^d
948 integer :: ipole,ic^d,inc^d,ineighbor,ipe_neighbor,ixs^
l,ixr^
l,n_i^d,idir
950 ipole=neighbor_pole(i^d,igrid)
953 {
do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
954 inc^db=2*i^db+ic^db\}
955 ipe_neighbor=neighbor_child(2,inc^d,igrid)
956 if(ipe_neighbor/=mype)
then
957 ineighbor=neighbor_child(1,inc^d,igrid)
964 shape=shape(psb(igrid)%ws(ixr^s,idir)))
970 {
do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
971 inc^db=2*i^db+ic^db\}
972 ipe_neighbor=neighbor_child(2,inc^d,igrid)
973 if(ipe_neighbor/=mype)
then
974 ineighbor=neighbor_child(1,inc^d,igrid)
977 n_i^d=i^d^d%n_i^dd=-i^dd;\}
987 shape=shape(psb(igrid)%ws(ixr^s,idir)))
988 call pole_copy_stg(psb(igrid)%ws,ixgs^ll,ixr^
l,pole_buf%ws,ixgs^ll,ixr^
l,idir,ipole)
995 end subroutine bc_fill_restrict_stg
998 subroutine bc_recv_prolong(igrid,i^D)
999 integer,
intent(in) :: igrid,i^d
1001 integer :: ic^d,ipe_neighbor,inc^d
1003 ic^d=1+modulo(node(pig^d_,igrid)-1,2);
1004 if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.})
return
1006 ipe_neighbor=neighbor(2,i^d,igrid)
1007 if (ipe_neighbor/=mype)
then
1010 itag=(3**^nd+4**^nd)*(igrid-1)+3**^nd+{inc^d*4**(^d-1)+}
1011 call mpi_irecv(psc(igrid)%w,1,
type_recv_p(inc^d), &
1013 if(stagger_grid)
then
1016 mpi_double_precision,ipe_neighbor,itag,&
1022 end subroutine bc_recv_prolong
1025 subroutine bc_send_prolong(igrid,i^D)
1026 integer,
intent(in) :: igrid,i^d
1028 integer :: ic^d,inc^d,n_i^d,n_inc^d,ineighbor,ipe_neighbor,ixs^
l,ipole,idir
1030 ipole=neighbor_pole(i^d,igrid)
1032 {
do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1033 inc^db=2*i^db+ic^db\}
1034 ipe_neighbor=neighbor_child(2,inc^d,igrid)
1035 if(ipe_neighbor/=mype)
then
1036 ineighbor=neighbor_child(1,inc^d,igrid)
1041 itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
1042 call mpi_isend(psb(igrid)%w,1,
type_send_p(inc^d), &
1044 if(stagger_grid)
then
1051 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
1052 ibuf_start=ibuf_next
1056 mpi_double_precision,ipe_neighbor,itag, &
1064 n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
1066 if(isend_buf(ipwbuf)/=0)
then
1069 deallocate(pwbuf(ipwbuf)%w)
1071 allocate(pwbuf(ipwbuf)%w(ixs^s,nwhead:nwtail))
1072 call pole_buffer(pwbuf(ipwbuf)%w,ixs^
l,ixs^
l,psb(igrid)%w,ixg^ll,ixs^
l,ipole)
1075 itag=(3**^nd+4**^nd)*(ineighbor-1)+3**^nd+{n_inc^d*4**(^d-1)+}
1076 isizes={(ixsmax^d-ixsmin^d+1)*}*nwbc
1077 call mpi_isend(pwbuf(ipwbuf)%w,isizes,mpi_double_precision, &
1079 ipwbuf=1+modulo(ipwbuf,
npwbuf)
1080 if(stagger_grid)
then
1087 reshape(psb(igrid)%ws(ixs^s,idir),shapes)
1088 ibuf_start=ibuf_next
1092 mpi_double_precision,ipe_neighbor,itag, &
1100 end subroutine bc_send_prolong
1103 subroutine bc_fill_prolong(igrid,i^D)
1104 integer,
intent(in) :: igrid,i^d
1106 integer :: ipe_neighbor,ineighbor,ixs^
l,ixr^
l,ic^d,inc^d,n_i^d,n_inc^d,ipole,idir
1108 ipole=neighbor_pole(i^d,igrid)
1111 {
do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1112 inc^db=2*i^db+ic^db\}
1113 ipe_neighbor=neighbor_child(2,inc^d,igrid)
1114 if(ipe_neighbor==mype)
then
1116 ineighbor=neighbor_child(1,inc^d,igrid)
1120 psc(ineighbor)%w(ixr^s,nwhead:nwtail) &
1121 =psb(igrid)%w(ixs^s,nwhead:nwtail)
1122 if(stagger_grid)
then
1126 psc(ineighbor)%ws(ixr^s,idir)=psb(igrid)%ws(ixs^s,idir)
1132 {
do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
1133 inc^db=2*i^db+ic^db\}
1134 ipe_neighbor=neighbor_child(2,inc^d,igrid)
1135 if(ipe_neighbor==mype)
then
1137 ineighbor=neighbor_child(1,inc^d,igrid)
1140 n_inc^d=inc^d^d%n_inc^dd=ic^dd-i^dd;\}
1143 call pole_copy(psc(ineighbor)%w,
ixcog^
l,ixr^
l,psb(igrid)%w,ixg^ll,ixs^
l,ipole)
1144 if(stagger_grid)
then
1148 call pole_copy_stg(psc(ineighbor)%ws,
ixcogs^
l,ixr^
l,psb(igrid)%ws,ixgs^ll,ixs^
l,idir,ipole)
1154 end subroutine bc_fill_prolong
1156 subroutine gc_prolong(igrid)
1157 integer,
intent(in) :: igrid
1159 integer :: i^d,idims,iside
1160 logical,
dimension(-1:1^D&) :: needprolong
1165 if (neighbor_type(i^d,igrid)==neighbor_coarse)
then
1166 call bc_prolong(igrid,i^d)
1167 needprolong(i^d)=.true.
1170 if(stagger_grid)
then
1181 if (needprolong(i^dd))
call bc_prolong_stg(igrid,i^dd,needprolong)
1192 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1198 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1204 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1210 if (needprolong(i^d))
call bc_prolong_stg(igrid,i^d,needprolong)
1213 end subroutine gc_prolong
1216 subroutine bc_fill_prolong_stg(igrid,i^D)
1217 integer,
intent(in) :: igrid,i^d
1219 integer :: ipe_neighbor,ineighbor,ipole,ixr^
l,ic^d,inc^d,n_inc^d,idir
1221 ic^d=1+modulo(node(pig^d_,igrid)-1,2);
1222 if ({.not.(i^d==0.or.i^d==2*ic^d-3)|.or.})
return
1224 ipe_neighbor=neighbor(2,i^d,igrid)
1225 if(ipe_neighbor/=mype)
then
1226 ineighbor=neighbor(1,i^d,igrid)
1227 ipole=neighbor_pole(i^d,igrid)
1236 shape=shape(psc(igrid)%ws(ixr^s,idir)))
1243 n_inc^d=2*i^d+(3-ic^d)^d%n_inc^dd=-2*i^dd+ic^dd;\}
1251 shape=shape(psc(igrid)%ws(ixr^s,idir)))
1252 call pole_copy_stg(psc(igrid)%ws,
ixcogs^
l,ixr^
l,pole_buf%ws,ixgs^ll,ixr^
l,idir,ipole)
1258 end subroutine bc_fill_prolong_stg
1261 subroutine bc_prolong(igrid,i^D)
1264 double precision :: dxfi^d, dxco^d, xfimin^d, xcomin^d, invdxco^d
1265 integer :: i^d,igrid
1266 integer :: ixfi^
l,ixco^
l,ii^d, idims,iside,ixb^
l
1269 dxfi^d=rnode(rpdx^d_,igrid);
1271 invdxco^d=1.d0/dxco^d;
1277 xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
1278 xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;
1280 if(phyboundblock(igrid).and.
bcphys)
then
1283 ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
1284 ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
1288 if(neighbor_type(-1,0,0,igrid)==neighbor_boundary .or. &
1289 neighbor_type(1,0,0,igrid)==neighbor_boundary)
then
1290 if(neighbor_type(0,-1,0,igrid)==neighbor_boundary) ixcomin2=ixcommin2
1291 if(neighbor_type(0,0,-1,igrid)==neighbor_boundary) ixcomin3=ixcommin3
1292 if(neighbor_type(0,1,0,igrid)==neighbor_boundary) ixcomax2=ixcommax2
1293 if(neighbor_type(0,0,1,igrid)==neighbor_boundary) ixcomax3=ixcommax3
1295 else if(idims == 2)
then
1296 if(neighbor_type(0,-1,0,igrid)==neighbor_boundary .or. &
1297 neighbor_type(0,1,0,igrid)==neighbor_boundary)
then
1298 if(neighbor_type(0,0,-1,igrid)==neighbor_boundary) ixcomin3=ixcommin3
1299 if(neighbor_type(0,0,1,igrid)==neighbor_boundary) ixcomax3=ixcommax3
1304 ii^d=kr(^d,idims)*(2*iside-3);
1305 if(neighbor_type(ii^d,igrid)/=neighbor_boundary) cycle
1306 if(( {(iside==1.and.idims==^d.and.ixcomin^d<ixcogmin^d+nghostcells)|.or.} ) &
1307 .or.( {(iside==2.and.idims==^d.and.ixcomax^d>ixcogmax^d-nghostcells)|.or. }))
then
1308 {ixbmin^d=merge(ixcogmin^d,ixcomin^d,idims==^d);}
1309 {ixbmax^d=merge(ixcogmax^d,ixcomax^d,idims==^d);}
1310 call bc_phys(iside,idims,time,0.d0,psc(igrid),
ixcog^
l,ixb^
l)
1316 if(prolongprimitive)
then
1323 ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
1324 ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
1329 call interpolation_copy(igrid,ixfi^
l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
1331 call interpolation_linear(igrid,ixfi^
l,dxfi^d,xfimin^d,dxco^d,invdxco^d,xcomin^d)
1334 if(prolongprimitive)
then
1339 end subroutine bc_prolong
1341 subroutine bc_prolong_stg(igrid,i^D,NeedProlong)
1343 double precision :: dxfi^d,dxco^d,xfimin^d,xcomin^d,invdxco^d
1344 integer :: igrid,i^d
1345 integer :: ixfi^
l,ixco^
l
1346 logical,
dimension(-1:1^D&) :: needprolong
1347 logical :: fine_^lin
1351 if(i^d>-1) fine_min^din=(.not.needprolong(i^dd-kr(^d,^dd)).and.neighbor_type(i^dd-kr(^d,^dd),igrid)/=1)
1352 if(i^d<1) fine_max^din=(.not.needprolong(i^dd+kr(^d,^dd)).and.neighbor_type(i^dd+kr(^d,^dd),igrid)/=1)
1357 dxfi^d=rnode(rpdx^d_,igrid);
1359 invdxco^d=1.d0/dxco^d;
1361 xfimin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxfi^d;
1362 xcomin^d=rnode(rpxmin^d_,igrid)-dble(nghostcells)*dxco^d;
1367 ixcomin^d=int((xfimin^d+(dble(ixfimin^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1-1;
1368 ixcomax^d=int((xfimin^d+(dble(ixfimax^d)-half)*dxfi^d-xcomin^d)*invdxco^d)+1+1;
1370 if(prolongprimitive)
call phys_to_primitive(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1372 call prolong_2nd_stg(psc(igrid),psb(igrid),ixco^
l,ixfi^
l,dxco^d,xcomin^d,dxfi^d,xfimin^d,.true.,fine_^lin)
1374 if(prolongprimitive)
call phys_to_conserved(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1377 needprolong(i^d)=.false.
1379 end subroutine bc_prolong_stg
1381 subroutine interpolation_linear(igrid,ixFi^L,dxFi^D,xFimin^D, &
1382 dxCo^D,invdxCo^D,xComin^D)
1384 integer,
intent(in) :: igrid, ixfi^
l
1385 double precision,
intent(in) :: dxfi^d, xfimin^d,dxco^d, invdxco^d, xcomin^d
1387 double precision :: xco^d, xfi^d, eta^d
1388 double precision :: slopel, sloper, slopec, signc, signr
1389 double precision :: slope(1:nw,ndim)
1391 double precision :: signedfactorhalf^d
1392 integer :: ixco^d, jxco^d, hxco^d, ixfi^d, ix^d, iw, idims, nwmin,nwmax
1397 if(prolongprimitive)
then
1405 {
do ixfi^db = ixfi^lim^db
1408 xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
1413 ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1
1417 xco^db=xcomin^db+(dble(ixco^db)-half)*dxco^db \}
1423 if(slab_uniform)
then
1432 eta^d=(xfi^d-xco^d)*invdxco^d;
1466 ix^d=2*int((ixfi^d+ixmlo^d)/2)-ixmlo^d;
1467 {
if(xfi^d>xco^d)
then
1468 signedfactorhalf^d=0.5d0
1470 signedfactorhalf^d=-0.5d0
1472 eta^d=signedfactorhalf^d*(one-psb(igrid)%dvolume(ixfi^dd) &
1473 /sum(psb(igrid)%dvolume(ix^d:ix^d+1^d%ixFi^dd))) \}
1480 hxco^d=ixco^d-kr(^d,idims)\
1481 jxco^d=ixco^d+kr(^d,idims)\
1484 slopel=psc(igrid)%w(ixco^d,iw)-psc(igrid)%w(hxco^d,iw)
1485 sloper=psc(igrid)%w(jxco^d,iw)-psc(igrid)%w(ixco^d,iw)
1486 slopec=half*(sloper+slopel)
1489 signr=sign(one,sloper)
1490 signc=sign(one,slopec)
1508 slope(iw,idims)=signc*max(zero,min(dabs(slopec), &
1509 signc*slopel,signc*sloper))
1515 psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)+&
1516 {(slope(nwmin:nwmax,^d)*eta^d)+}
1520 if(prolongprimitive)
then
1522 call phys_to_conserved(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1525 end subroutine interpolation_linear
1527 subroutine interpolation_copy(igrid, ixFi^L,dxFi^D,xFimin^D, &
1528 dxCo^D,invdxCo^D,xComin^D)
1530 integer,
intent(in) :: igrid, ixfi^
l
1531 double precision,
intent(in) :: dxfi^d, xfimin^d,dxco^d, invdxco^d, xcomin^d
1533 double precision :: xfi^d
1534 integer :: ixco^d, ixfi^d, nwmin,nwmax
1536 if(prolongprimitive)
then
1544 {
do ixfi^db = ixfi^lim^db
1546 xfi^db=xfimin^db+(dble(ixfi^db)-half)*dxfi^db
1550 ixco^db=int((xfi^db-xcomin^db)*invdxco^db)+1\}
1553 psb(igrid)%w(ixfi^d,nwmin:nwmax)=psc(igrid)%w(ixco^d,nwmin:nwmax)
1557 if(prolongprimitive)
call phys_to_conserved(ixg^ll,ixfi^
l,psb(igrid)%w,psb(igrid)%x)
1559 end subroutine interpolation_copy
1561 subroutine pole_copy(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,ipole)
1563 integer,
intent(in) :: ixir^
l,ixr^
l,ixis^
l,ixs^
l,ipole
1564 double precision :: wrecv(ixir^s,1:nw), wsend(ixis^s,1:nw)
1566 integer :: iw, iside, ib
1570 iside=int((i^d+3)/2)
1573 select case (typeboundary(iw,ib))
1575 wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1577 wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1579 call mpistop(
"Pole boundary condition should be symm or asymm")
1584 end subroutine pole_copy
1586 subroutine pole_copy_stg(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,idirs,ipole)
1588 integer,
intent(in) :: ixir^
l,ixr^
l,ixis^
l,ixs^
l,idirs,ipole
1590 double precision :: wrecv(ixir^s,1:nws), wsend(ixis^s,1:nws)
1591 integer :: ib, iside
1595 iside=int((i^d+3)/2)
1597 select case (typeboundary(iw_mag(idirs),ib))
1599 wrecv(ixr^s,idirs) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
1601 wrecv(ixr^s,idirs) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,idirs)
1603 call mpistop(
"Pole boundary condition should be symm or asymm")
1608 end subroutine pole_copy_stg
1610 subroutine pole_buffer(wrecv,ixIR^L,ixR^L,wsend,ixIS^L,ixS^L,ipole)
1612 integer,
intent(in) :: ixir^
l,ixr^
l,ixis^
l,ixs^
l,ipole
1613 double precision :: wrecv(ixir^s,nwhead:nwtail), wsend(ixis^s,1:nw)
1615 integer :: iw, iside, ib
1619 iside=int((i^d+3)/2)
1622 select case (typeboundary(iw,ib))
1624 wrecv(ixr^s,iw) = wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1626 wrecv(ixr^s,iw) =-wsend(ixsmax^d:ixsmin^d:-1^d%ixS^s,iw)
1628 call mpistop(
"Pole boundary condition should be symm or asymm")
1633 end subroutine pole_buffer