boundary_conditions.t

Go to the documentation of this file.

!> fill ghost cells at a physical boundary
subroutine bc_phys(iside,idims,time,qdt,s,ixG^L,ixB^L)
  use mod_usr_methods, only: usr_special_bc
  use mod_bc_data, only: bc_data_set
  use mod_global_parameters
  use mod_physics
  use mod_mhd_phys, only: get_divb
 
  integer, intent(in) :: iside, idims, ixG^L,ixB^L
  double precision, intent(in) :: time,qdt
  type(state), intent(inout) :: s
  double precision :: wtmp(ixG^S,1:nwflux)
 
  integer :: idir, is
  integer :: ixOs^L,hxO^L,jxO^L
  double precision :: Q(ixG^S),Qp(ixG^S) 
  integer :: iw, iB, ix^D, ixO^L, ixM^L, nghostcellsi,iib^D
  logical  :: isphysbound
 
  associate(x=>s%x,w=>s%w,ws=>s%ws)
  select case (idims)
  {case (^d)
     if (iside==2) then
        ! maximal boundary
        ib=2*^d
        ixomin^dd=ixbmax^d+1-nghostcells^d%ixOmin^dd=ixbmin^dd;
        ixomax^dd=ixbmax^dd;
        ! cont/symm/asymm types
        do iw=1,nwflux+nwaux
           select case (typeboundary(iw,ib))
           case (bc_special)
              ! skip it here, do AFTER all normal type boundaries are set
           case (bc_cont)
              do ix^d=ixomin^d,ixomax^d
                 w(ix^d^d%ixO^s,iw) = w(ixomin^d-1^d%ixO^s,iw)
              end do
           case (bc_symm)
              w(ixo^s,iw) = w(ixomin^d-1:ixomin^d-nghostcells:-1^d%ixO^s,iw)
           case (bc_asymm)
              w(ixo^s,iw) =-w(ixomin^d-1:ixomin^d-nghostcells:-1^d%ixO^s,iw)
           case (bc_periodic)
              ! skip it here, periodic bc info should come from neighbors
           case(bc_noinflow)
              if (iw==1+^d)then
                do ix^d=ixomin^d,ixomax^d
                    w(ix^d^d%ixO^s,iw) = max(w(ixomin^d-1^d%ixO^s,iw),zero)
                end do
              else
                do ix^d=ixomin^d,ixomax^d
                    w(ix^d^d%ixO^s,iw) = w(ixomin^d-1^d%ixO^s,iw)
                end do
              end if
           case (bc_aperiodic)
              !this just multiplies the variables with (-), they have been set from neighbors just like periodic.
              w(ixo^s,iw) = - w(ixo^s,iw)
           case (bc_data)
              ! skip it here, do AFTER all normal type boundaries are set
           case (bc_character)
              ! skip it here, do AFTER all normal type boundaries are set
           case default
              write (unitterm,*) "Undefined boundarytype found in bc_phys", &
                 "for variable iw=",iw," and side iB=",ib
           end select
        end do
        if(stagger_grid) then
          do idir=1,nws
          ! At this stage, extrapolation is applied only to the tangential components
            if(idir==^d) cycle 
            ixosmax^dd=ixomax^dd;
            ixosmin^dd=ixomin^dd-kr(^dd,idir);
            select case(typeboundary(iw_mag(idir),ib))
            case (bc_special)
               ! skip it here, do AFTER all normal type boundaries are set
            case (bc_cont)
              do ix^d=ixosmin^d,ixosmax^d
                 ws(ix^d^d%ixOs^s,idir) = ws(ixosmin^d-1^d%ixOs^s,idir)
              end do
            case (bc_symm)
              ws(ixos^s,idir) = ws(ixosmin^d-1:ixosmin^d-nghostcells:-1^d%ixOs^s,idir)
            case (bc_asymm)
              ws(ixos^s,idir) =-ws(ixosmin^d-1:ixosmin^d-nghostcells:-1^d%ixOs^s,idir)
            case (bc_periodic)
            case default
              write (unitterm,*) "Undefined boundarytype found in bc_phys", &
                 "for variable iw=",iw," and side iB=",ib
            end select
          end do
          ! Now that the tangential components are set,
          ! fill the normal components using a prescription for the divergence.
          ! This prescription is given by the typeB for the normal component.
          do idir=1,nws
            ! Consider only normal direction
            if (idir/=^d) cycle
            ixos^l=ixo^l;
            select case(typeboundary(iw_mag(idir),ib))
            case(bc_cont)
              hxo^l=ixo^l-nghostcells*kr(^dd,^d);
              ! Calculate divergence and partial divergence
              call get_divb(s%w,ixg^l,hxo^l,q)
              ws(ixos^s,idir)=zero
              do ix^d=0,nghostcells-1
                call get_divb(s%w,ixg^l,ixo^l,qp)
                ws(ixosmin^d+ix^d^d%ixOs^s,idir)=&
                  (q(hxomax^d^d%hxO^s)*s%dvolume(hxomax^d^d%hxO^s)&
                 -qp(ixomin^d+ix^d^d%ixO^s)*s%dvolume(ixomin^d+ix^d^d%ixO^s))&
                  /s%surfaceC(ixosmin^d+ix^d^d%ixOs^s,^d)
              end do
            case(bc_symm)
              ! (a)symmetric normal B ensures symmetric divB
              ws(ixos^s,idir)= ws(ixosmin^d-2:ixosmin^d-nghostcells-1:-1^d%ixOs^s,idir)
            case(bc_asymm)
              ! (a)symmetric normal B ensures symmetric divB
              ws(ixos^s,idir)=-ws(ixosmin^d-2:ixosmin^d-nghostcells-1:-1^d%ixOs^s,idir)
            case(bc_periodic)
            end select
          end do
          ! Fill cell averages
          call phys_face_to_center(ixo^l,s)
        end if
     else
        ! minimal boundary
        ib=2*^d-1
        ixomin^dd=ixbmin^dd;
        ixomax^dd=ixbmin^d-1+nghostcells^d%ixOmax^dd=ixbmax^dd;
        ! cont/symm/asymm types
        do iw=1,nwflux+nwaux
           select case (typeboundary(iw,ib))
           case (bc_special)
              ! skip it here, do AFTER all normal type boundaries are set
           case (bc_cont)
              do ix^d=ixomin^d,ixomax^d
                 w(ix^d^d%ixO^s,iw) = w(ixomax^d+1^d%ixO^s,iw)
              end do
           case (bc_symm)
              w(ixo^s,iw) = w(ixomax^d+nghostcells:ixomax^d+1:-1^d%ixO^s,iw)
           case (bc_asymm)
              w(ixo^s,iw) =-w(ixomax^d+nghostcells:ixomax^d+1:-1^d%ixO^s,iw)
           case (bc_periodic)
              ! skip it here, periodic bc info should come from neighbors
           case(bc_noinflow)
              if (iw==1+^d)then
                 do ix^d=ixomin^d,ixomax^d
                   w(ix^d^d%ixO^s,iw) = min(w(ixomax^d+1^d%ixO^s,iw),zero)
                 end do
              else
                 do ix^d=ixomin^d,ixomax^d
                   w(ix^d^d%ixO^s,iw) = w(ixomax^d+1^d%ixO^s,iw)
                 end do
              end if
           case (bc_aperiodic)
              !this just multiplies the variables with (-), they have been set from neighbors just like periodic.
              w(ixo^s,iw) = - w(ixo^s,iw)
           case (bc_data)
              ! skip it here, do AFTER all normal type boundaries are set
           case (bc_character)
              ! skip it here, do AFTER all normal type boundaries are set
           case default
              write (unitterm,*) "Undefined boundarytype found in bc_phys", &
                 "for variable iw=",iw," and side iB=",ib
           end select
        end do
        if(stagger_grid) then
          do idir=1,nws
          ! At this stage, extrapolation is applied only to the tangential components
            if(idir==^d) cycle 
            ixosmax^dd=ixomax^dd;
            ixosmin^dd=ixomin^dd-kr(^dd,idir);
            select case(typeboundary(iw_mag(idir),ib))
            case (bc_special)
              ! skip it here, periodic bc info should come from neighbors
            case (bc_cont)
              do ix^d=ixosmin^d,ixosmax^d
                 ws(ix^d^d%ixOs^s,idir) = ws(ixosmax^d+1^d%ixOs^s,idir)
              end do
            case (bc_symm)
              ws(ixos^s,idir) = ws(ixosmax^d+nghostcells:ixosmax^d+1:-1^d%ixOs^s,idir)
            case (bc_asymm)
              ws(ixos^s,idir) =-ws(ixosmax^d+nghostcells:ixosmax^d+1:-1^d%ixOs^s,idir)
            case (bc_periodic)
            case default
              write (unitterm,*) "Undefined boundarytype in bc_phys", &
                 "for variable iw=",iw," and side iB=",ib
            end select
          end do
          ! Now that the tangential components are set,
          ! fill the normal components using a prescription for the divergence.
          ! This prescription is given by the typeB for the normal component.
          do idir=1,nws
            ! Consider only normal direction
            if (idir/=^d) cycle
            ixos^l=ixo^l-kr(^dd,^d);
            select case(typeboundary(iw_mag(idir),ib))
            case(bc_cont)
              jxo^l=ixo^l+nghostcells*kr(^dd,^d);
              ! Calculate divergence and partial divergence
              call get_divb(s%w,ixg^l,jxo^l,q)
              ws(ixos^s,idir)=zero
              do ix^d=0,nghostcells-1
                call get_divb(s%w,ixg^l,ixo^l,qp)
                ws(ixosmax^d-ix^d^d%ixOs^s,idir)=&
                 -(q(jxomin^d^d%jxO^s)*s%dvolume(jxomin^d^d%jxO^s)&
                 -qp(ixomax^d-ix^d^d%ixO^s)*s%dvolume(ixomax^d-ix^d^d%ixO^s))&
                  /s%surfaceC(ixosmax^d-ix^d^d%ixOs^s,^d)
              end do
            case(bc_symm)
              ! (a)symmetric normal B ensures symmetric divB
              ws(ixos^s,idir)= ws(ixosmax^d+nghostcells+1:ixosmax^d+2:-1^d%ixOs^s,idir)
            case(bc_asymm)
              ! (a)symmetric normal B ensures symmetric divB
              ws(ixos^s,idir)=-ws(ixosmax^d+nghostcells+1:ixosmax^d+2:-1^d%ixOs^s,idir)
            case(bc_periodic)
            end select
          end do
          ! Fill cell averages
          call phys_face_to_center(ixo^l,s)
        end if
     end if \}
  end select
 
  ! do user defined special boundary conditions
  if (any(typeboundary(1:nwflux+nwaux,ib)==bc_special)) then
     if (.not. associated(usr_special_bc)) &
          call mpistop("usr_special_bc not defined")
     call usr_special_bc(time,ixg^l,ixo^l,ib,w,x)
  end if
 
  ! fill boundary conditions from external data vtk files
  if (any(typeboundary(1:nwflux+nwaux,ib)==bc_data)) then
     call bc_data_set(time,ixg^l,ixo^l,ib,w,x)
  end if
 
  {#IFDEF EVOLVINGBOUNDARY
  if (any(typeboundary(1:nwflux,ib)==bc_character)) then
    ixm^l=ixm^ll;
    if(ixgmax1==ixghi1) then
      nghostcellsi=nghostcells
    else
      nghostcellsi=ceiling(nghostcells*0.5d0)
    end if
    select case (idims)
    {case (^d)
       if (iside==2) then
          ! maximal boundary
          ixomin^dd=ixgmax^d+1-nghostcellsi^d%ixOmin^dd=ixbmin^dd;
          ixomax^dd=ixbmax^dd;
          if(all(w(ixo^s,1:nwflux)==0.d0)) then
            do ix^d=ixomin^d,ixomax^d
               w(ix^d^d%ixO^s,1:nwflux) = w(ixomin^d-1^d%ixO^s,1:nwflux)
            end do
          end if
          if(qdt>0.d0.and.ixgmax^d==ixghi^d) then
            ixomin^dd=ixomin^d^d%ixOmin^dd=ixmmin^dd;
            ixomax^dd=ixomax^d^d%ixOmax^dd=ixmmax^dd;
            wtmp(ixg^s,1:nw)=pso(block%igrid)%w(ixg^s,1:nw)
            call characteristic_project(idims,iside,ixg^l,ixo^l,wtmp,x,dxlevel,qdt)
            w(ixo^s,1:nwflux)=wtmp(ixo^s,1:nwflux)
          end if
       else
          ! minimal boundary
          ixomin^dd=ixbmin^dd;
          ixomax^dd=ixgmin^d-1+nghostcellsi^d%ixOmax^dd=ixbmax^dd;
          if(all(w(ixo^s,1:nwflux)==0.d0)) then
            do ix^d=ixomin^d,ixomax^d
               w(ix^d^d%ixO^s,1:nwflux) = w(ixomax^d+1^d%ixO^s,1:nwflux)
            end do
          end if
          if(qdt>0.d0.and.ixgmax^d==ixghi^d) then
            ixomin^dd=ixomin^d^d%ixOmin^dd=ixmmin^dd;
            ixomax^dd=ixomax^d^d%ixOmax^dd=ixmmax^dd;
            wtmp(ixg^s,1:nw)=pso(block%igrid)%w(ixg^s,1:nw)
            call characteristic_project(idims,iside,ixg^l,ixo^l,wtmp,x,dxlevel,qdt)
            w(ixo^s,1:nwflux)=wtmp(ixo^s,1:nwflux)
          end if
       end if \}
    end select
    if(ixgmax1==ixghi1) then
      call identifyphysbound(block%igrid,isphysbound,iib^d)   
      if(iib1==-1.and.iib2==-1) then
        do ix2=nghostcells,1,-1 
          do ix1=nghostcells,1,-1 
            w(ix^d,1:nwflux)=(w(ix1+1,ix2+1,1:nwflux)+w(ix1+1,ix2,1:nwflux)+w(ix1,ix2+1,1:nwflux))/3.d0
          end do
        end do
      end if
      if(iib1== 1.and.iib2==-1) then
        do ix2=nghostcells,1,-1 
          do ix1=ixmmax1+1,ixgmax1
            w(ix^d,1:nwflux)=(w(ix1-1,ix2+1,1:nwflux)+w(ix1-1,ix2,1:nwflux)+w(ix1,ix2+1,1:nwflux))/3.d0
          end do
        end do
      end if
      if(iib1==-1.and.iib2== 1) then
        do ix2=ixmmax2+1,ixgmax2
          do ix1=nghostcells,1,-1 
            w(ix^d,1:nwflux)=(w(ix1+1,ix2-1,1:nwflux)+w(ix1+1,ix2,1:nwflux)+w(ix1,ix2-1,1:nwflux))/3.d0
          end do
        end do
      end if
      if(iib1== 1.and.iib2== 1) then
        do ix2=ixmmax2+1,ixgmax2
          do ix1=ixmmax1+1,ixgmax1
            w(ix^d,1:nwflux)=(w(ix1-1,ix2-1,1:nwflux)+w(ix1-1,ix2,1:nwflux)+w(ix1,ix2-1,1:nwflux))/3.d0
          end do
        end do
      end if
    end if
  end if
  }
  !end do
 
  end associate
end subroutine bc_phys
 
!> fill inner boundary values
subroutine getintbc(time,ixG^L)
  use mod_usr_methods, only: usr_internal_bc
  use mod_global_parameters
 
  double precision, intent(in)   :: time
  integer, intent(in)            :: ixG^L
 
  integer :: iigrid, igrid, ixO^L
 
  ixo^l=ixg^l^lsubnghostcells;
 
  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
     block=>ps(igrid)
     ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
 
     if (associated(usr_internal_bc)) then
        call usr_internal_bc(node(plevel_,igrid),time,ixg^l,ixo^l,ps(igrid)%w,ps(igrid)%x)
     end if
  end do
  !$OMP END PARALLEL DO
 
end subroutine getintbc