mod_boundary_conditions.t

Go to the documentation of this file.

module mod_boundary_conditions
 
  implicit none
  private
 
  public :: bc_phys
  public :: getintbc
 
contains
 
  !> 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_functions_bfield, only: get_divb
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: iside, idims, ixg^l,ixb^l
    double precision, intent(in) :: time,qdt
    type(state), intent(inout) :: s
 
    double precision :: q(ixg^s),qp(ixg^s) 
    integer :: idir, ixos^l,hxo^l,jxo^l
    integer :: iw, ib, ix^d, ixo^l
 
    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_icarus)
              ! 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
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              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)
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              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)
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              case(bc_periodic)
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              end select
            end do
          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_icarus)
              ! 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
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              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)
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              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)
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              case(bc_periodic)
                ! Fill cell averages
                call phys_face_to_center(ixo^l,s)
              end select
            end do
          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
 
    ! fill boundary conditions from external data vtk files and do user defined special boundary conditions
    if (any(typeboundary(1:nwflux+nwaux,ib)==bc_icarus)) then
       call bc_data_set(time,ixg^l,ixo^l,ib,w,x)
       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
 
    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
 
end module mod_boundary_conditions