mod_point_searching.t

Go to the documentation of this file.

! for Cartesian coordinate system
module mod_point_searching
  use mod_global_parameters
  use mod_physics
  use mod_particle_base
  use mod_comm_lib, only: mpistop
  implicit none
 
contains
 
  subroutine get_point_w_ingrid(igrid,xp,wp,variable_type)
    double precision :: xp(1:ndim),wp(1:nw)
    character(*) :: variable_type
    integer, intent(in) :: igrid
 
    double precision :: dxb^D,xd^D,xb^L,temp
    integer :: ixI^L,ixO^L
    integer :: ixbl^D,ix^D,ixA^L,j,ingrid
    double precision :: factor(0:1^D&)
 
    ixi^l=ixg^ll;
    ixo^l=ixm^ll;
    ^d&xbmin^d=rnode(rpxmin^d_,igrid)-rnode(rpdx^d_,igrid);
    ^d&xbmax^d=rnode(rpxmax^d_,igrid)+rnode(rpdx^d_,igrid);
    ingrid=0
    {if (xp(^db)>xbmin^db .and. xp(^db)<xbmax^db) ingrid=ingrid+1\}
 
    if (ingrid==ndim) then
      ^d&dxb^d=rnode(rpdx^d_,igrid)\
      ^d&ixbl^d=floor((xp(^d)-ps(igrid)%x(ixomin^dd,^d))/dxb^d)+ixomin^d\
      ^d&xd^d=(xp(^d)-ps(igrid)%x(ixbl^dd,^d))/dxb^d\
      ^d&ixamin^d=ixbl^d\
      ^d&ixamax^d=ixbl^d+1\
 
      {do ix^d=0,1\}
        factor(ix^d)={abs(1-ix^d-xd^d)*}
      {enddo\}
 
      select case(variable_type)
        case('primitive')
          call phys_to_primitive(ixi^l,ixa^l,ps(igrid)%w,ps(igrid)%x)
        case('conserved')
 
        case default
          call mpistop("get_point_w_ingrid: Unknown variable type!")
      end select
 
      wp=0.d0
      {do ix^d=0,1\}
        do j=1,nw
          wp(j)=wp(j)+factor(ix^d)*ps(igrid)%w(ixbl^d+ix^d,j)
        enddo
      {enddo\}
 
      select case(variable_type)
        case('primitive')
        call phys_to_conserved(ixi^l,ixa^l,ps(igrid)%w,ps(igrid)%x)
        case('conserved')
 
        case default
          call mpistop("get_point_w_ingrid: Unknown variable type!")
      end select
 
      if(physics_type=='mhd') then
        wp(iw_mag(1):iw_mag(ndir))=0.d0
        {do ix^d=0,1\}
          do j=1,ndir
            if (b0field) then
              temp=ps(igrid)%w(ixbl^d+ix^d,iw_mag(j))+&
                                    ps(igrid)%B0(ixbl^d+ix^d,j,0)
              wp(iw_mag(j))=wp(iw_mag(j))+factor(ix^d)*temp
            else
              temp=ps(igrid)%w(ixbl^d+ix^d,iw_mag(j))
              wp(iw_mag(j))=wp(iw_mag(j))+factor(ix^d)*temp
            endif
          enddo
        {enddo\}
      endif
    else
      call mpistop("get_point_w_ingrid: The point is not in given grid!")
    endif
 
  end subroutine get_point_w_ingrid
 
  subroutine get_point_w(xp,wp,variable_type)
    ! for given point (xp), provide the plasma parameters (wp) at this point
 
    double precision :: xp(1:ndim),wp(1:nw)
    character(*) :: variable_type
 
    double precision :: x3d(3)
    integer :: indomain,ipe,igrid,j
 
    indomain=0
    {if (xp(^db)>=xprobmin^db .and. xp(^db)<xprobmax^db) indomain=indomain+1\}
    if (indomain==ndim) then
 
      ! find pe and igrid
      x3d=0.d0
      do j=1,ndim
        x3d(j)=xp(j)
      enddo
      call find_particle_ipe(x3d,igrid,ipe)
 
      !do interpolation to get point values
      if (mype==ipe) then
        call get_point_w_ingrid(igrid,xp,wp,variable_type)
      endif
 
      call mpi_bcast(wp,nw,mpi_double_precision,ipe,icomm,ierrmpi)
    endif
 
  end subroutine get_point_w
 
  subroutine get_cell_w(xp,wc,variable_type)
    ! for given point (xp), looking for corresponding cell and then provide
    ! the plasma parameters (wc) of this cell
 
    double precision :: xp(1:ndim),wc(1:nw)
    character(*) :: variable_type
 
    double precision :: x3d(3)
    double precision :: dxb^D,xb^L
    integer :: indomain,ixO^L,ixb^D,ixA^L,ixI^L
    integer :: ipe,igrid,j
 
    indomain=0
    {if (xp(^db)>=xprobmin^db .and. xp(^db)<xprobmax^db) indomain=indomain+1\}
    if (indomain==ndim) then
 
      ! find pe and igrid
      x3d=0.d0
      do j=1,ndim
        x3d(j)=xp(j)
      enddo
      call find_particle_ipe(x3d,igrid,ipe)
 
      if (mype==ipe) then
        ! looking for the cell index
        ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
        ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
        ^d&dxb^d=rnode(rpdx^d_,igrid)\
        ^d&ixomin^d=ixmlo^d\
        ^d&iximin^d=ixglo^d\
        ^d&iximax^d=ixghi^d\
        ^d&ixb^d=floor((xp(^d)-xbmin^d)/dxb^d)+ixomin^d\
        ^d&ixamin^d=ixb^d\
        ^d&ixamax^d=ixb^d\
 
        wc=0.d0
        select case(variable_type)
          case('primitive')
            call phys_to_primitive(ixi^l,ixa^l,ps(igrid)%w,ps(igrid)%x)
            do j=1,nw
              wc(j)=ps(igrid)%w(ixb^d,j)
            enddo
            call phys_to_conserved(ixi^l,ixa^l,ps(igrid)%w,ps(igrid)%x)
          case('conserved')
            do j=1,nw
              wc(j)=ps(igrid)%w(ixb^d,j)
            enddo
          case default
            call mpistop("get_point_w: Unknown variable type!")
        end select
 
        ! for magnetic field
        if(physics_type=='mhd') then
          wc(iw_mag(1):iw_mag(ndir))=0.d0
          do j=1,ndir
            if (b0field) then
              wc(iw_mag(j))=ps(igrid)%w(ixb^d,iw_mag(j))+&
                            ps(igrid)%B0(ixb^d,j,0)
            else
              wc(iw_mag(j))=ps(igrid)%w(ixb^d,iw_mag(j))
            endif
          enddo
        endif
      endif
 
      call mpi_bcast(wc,nw,mpi_double_precision,ipe,icomm,ierrmpi)
    endif
 
  end subroutine get_cell_w
 
  subroutine get_cell_index(xp,ipe,igrid,ixc^D)
    ! for given point (xp), provide the corrsponding processor (ipe), 
    ! grid number (igrid) and cell index (ixc^D)
 
    double precision :: xp(1:ndim)
    integer :: ipe,igrid,ixc^D
 
    double precision :: x3d(1:3)
    double precision :: dxb^D,xb^L
    integer :: indomain,ixO^L,j
    integer :: datas(1:ndim+2)
 
    indomain=0
    {if (xp(^db)>=xprobmin^db .and. xp(^db)<xprobmax^db) indomain=indomain+1\}
 
    if (indomain==ndim) then
      ! find pe and igrid
      x3d=0.d0
      do j=1,ndim
        x3d(j)=xp(j)
      enddo
      call find_particle_ipe(x3d,igrid,ipe)
 
      if (mype==ipe) then
        ! looking for the cell index
        ^d&xbmin^d=rnode(rpxmin^d_,igrid)\
        ^d&xbmax^d=rnode(rpxmax^d_,igrid)\
        ^d&dxb^d=rnode(rpdx^d_,igrid)\
        ^d&ixomin^d=ixmlo^d\
        ^d&ixc^d=floor((xp(^d)-xbmin^d)/dxb^d)+ixomin^d\
        ^d&datas(^d)=ixc^d\
        datas(ndim+1)=ipe
        datas(ndim+2)=igrid
      endif
 
      call mpi_bcast(datas,ndim+2,mpi_integer,ipe,icomm,ierrmpi)
    endif
 
    ^d&ixc^d=datas(^d)\
    ipe=datas(ndim+1)
    igrid=datas(ndim+2)
 
  end subroutine get_cell_index
 
end module mod_point_searching