get_igslice Subroutine

public subroutine get_igslice(dir, xphys, igslice)

Uses

  • proc~~get_igslice~~UsesGraph proc~get_igslice get_igslice module~mod_geometry mod_geometry proc~get_igslice->module~mod_geometry module~mod_global_parameters mod_global_parameters proc~get_igslice->module~mod_global_parameters module~mod_comm_lib mod_comm_lib module~mod_geometry->module~mod_comm_lib module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_variables->module~mod_basic_types

Block indices, per level, of the blocks a slice passes through.

The position comes in as a physical coordinate, because that is what the rest of the slice machinery uses - fill_subnode picks the cell within a block by comparing against ps(igrid)%x. Everything below, though, is block-index arithmetic in the logical coordinate, the one xprobmin, dg and ng are expressed in. Under LOG_RADIUS those differ, so the conversion is done here, once, and the body is unchanged.

Getting this backwards is a silent trap: converting the slice position at read time instead would fix this routine and break fill_subnode, whose minloc would then match nothing, return 0, and emit a slice with all-zero coordinates without any error.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: dir
double precision, intent(in) :: xphys

slice position, a physical coordinate
integer, intent(out), dimension(nlevelshi) :: igslice

Calls

proc~~get_igslice~~CallsGraph proc~get_igslice get_igslice proc~mpistop mpistop proc~get_igslice->proc~mpistop mpi_abort mpi_abort proc~mpistop->mpi_abort

Called by

proc~~get_igslice~~CalledByGraph proc~get_igslice get_igslice proc~find_particle_ipe find_particle_ipe proc~find_particle_ipe->proc~get_igslice proc~select_slice select_slice proc~select_slice->proc~get_igslice proc~advect_create_particles advect_create_particles proc~advect_create_particles->proc~find_particle_ipe proc~apply_periodb apply_periodB proc~apply_periodb->proc~find_particle_ipe proc~comm_particles comm_particles proc~comm_particles->proc~find_particle_ipe proc~comm_particles->proc~apply_periodb proc~comm_particles_global comm_particles_global proc~comm_particles_global->proc~find_particle_ipe proc~gca_create_particles gca_create_particles proc~gca_create_particles->proc~find_particle_ipe proc~get_cell_index get_cell_index proc~get_cell_index->proc~find_particle_ipe proc~get_cell_w get_cell_w proc~get_cell_w->proc~find_particle_ipe proc~get_point_w get_point_w proc~get_point_w->proc~find_particle_ipe proc~lorentz_create_particles Lorentz_create_particles proc~lorentz_create_particles->proc~find_particle_ipe proc~put_slice put_slice proc~put_slice->proc~select_slice proc~read_from_snapshot read_from_snapshot proc~read_from_snapshot->proc~find_particle_ipe proc~sample_create_particles sample_create_particles proc~sample_create_particles->proc~find_particle_ipe proc~advance_particles advance_particles proc~advance_particles->proc~comm_particles proc~handle_particles handle_particles proc~handle_particles->proc~comm_particles_global proc~handle_particles->proc~advance_particles proc~particles_create particles_create proc~particles_create->proc~advect_create_particles proc~particles_create->proc~gca_create_particles proc~particles_create->proc~lorentz_create_particles proc~particles_create->proc~sample_create_particles proc~read_particles_snapshot read_particles_snapshot proc~read_particles_snapshot->proc~comm_particles_global proc~read_particles_snapshot->proc~read_from_snapshot proc~write_slice write_slice proc~write_slice->proc~put_slice proc~advance advance proc~advance->proc~handle_particles proc~main main proc~main->proc~handle_particles proc~main->proc~particles_create proc~main->proc~read_particles_snapshot proc~timeintegration timeintegration proc~main->proc~timeintegration proc~saveamrfile saveamrfile proc~saveamrfile->proc~write_slice proc~timeintegration->proc~advance proc~timeintegration->proc~saveamrfile program~agile agile program~agile->proc~main

Source Code

  subroutine get_igslice(dir,xphys,igslice)
    use mod_global_parameters
    use mod_geometry, only: s_of_r
    integer, intent(in) :: dir
    !> slice position, a physical coordinate
    double precision, intent(in) :: xphys
    integer, dimension(nlevelshi), intent(out) :: igslice
    ! .. local ..
    integer :: level
    double precision :: distance
    !> the same position in the logical coordinate used below
    double precision :: x
    !double precision :: xsgrid(ndim,nlevelshi),qs(ndim),xmgrid(ndim,3),xnew,xlgrid(ndim,3)
    !integer :: nbefore

    if (xphys.ne.xphys) call mpistop("get_igslice: your slice position is NaN!")

    x = xphys

    select case (dir)
    case (1)
      select case (stretch_type(1))
      case (stretch_none) ! This dimension is unstretched
        if (x<=xprobmin1) then
          igslice=1
        else if (x>=xprobmax1) then
          do level=1,refine_max_level
            igslice(level)=ng1(level)
          end do
        else
          distance=x-xprobmin1
          do level = 1, refine_max_level
            igslice(level) = ceiling(distance/dg1(level))
          end do
        end if

      case (stretch_uni) ! Uniform stretching

        if (x<=xprobmin1) then
          igslice=1
        else if (x>=xprobmax1) then
          do level=1,refine_max_level
            igslice(level)=ng1(level)
          end do
        else
          distance=x-xprobmin1
          do level=1,refine_max_level
            igslice(level)= ceiling(dlog(distance/dxfirst(level,&
               1)*(qstretch(level,1)-1.d0)+1.d0)/dlog(qstretch(level,&
               1))/dble(block_nx1))
          end do
        end if

      case (stretch_symm) ! Symmetric stretching
        ! symmetric stretch about 0.5*(xprobmin+xprobmax)
        if (x<=xprobmin1) then
          igslice=1
        else if (x>=xprobmax1) then
          do level=1,refine_max_level
            igslice(level)=ng1(level)
          end do
        else if(x<xprobmin1+xstretch1) then
          distance=xprobmin1+xstretch1-x
          ! stretch to left from xprobmin+xstretch
          do level=1,refine_max_level
            igslice(level) = nstretchedblocks(level,&
               1)/2-int(dlog(distance*(qstretch(level,1)-one)/dxfirst(level,&
               1)+one)/dlog(qstretch(level,1))/dble(block_nx1))
          end do
        else if(x>xprobmax1-xstretch1) then
          distance=x-xprobmax1+xstretch1
          ! stretch to right from xprobmax-xstretch
          do level=1,refine_max_level
            igslice(level) = ceiling(dlog(distance*(qstretch(level,&
               1)-one)/dxfirst(level,1)+one)/dlog(qstretch(level,&
               1))/dble(block_nx1))+ng1(level)-nstretchedblocks(level,1)/2
          end do
        else
          ! possible non-stretched central part
          distance=x-xprobmin1-xstretch1
          do level=1,refine_max_level
            igslice(level)=nstretchedblocks(level,&
               1)/2+ceiling(distance/dg1(level))
          end do
        end if
 !xsgrid(1,1)=half*(xprobmax1-xprobmin1)!        /(half*domain_nx1-half*nstretchedblocks_baselevel(1)*block_nx1 !        +(one-qstretch_baselevel(1)**(half*nstretchedblocks_baselevel(1)*block_nx1)) !        /(one-qstretch_baselevel(1)))
 !xmgrid(1,1)=xsgrid(1,1)*(domain_nx1-nstretchedblocks_baselevel(1)*block_nx1)
        !xlgrid(1,1)=xsgrid(1,1)

        !if (x .ge. xmgrid(1,1)*half) then
        !  xnew=x-xmgrid(1,1)*half
        !  do level=1,refine_max_level
        !    qs(1)=qstretch_baselevel(1)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(1,level)=xsgrid(1,level-1)/(one+qs(1))
 !nbefore=(domain_nx1/block_nx1-nstretchedblocks_baselevel(1)/2)*2**(level-1)
 !igslice(level)=floor((dlog(one-xnew/xsgrid(1,level)*(one-qs(1)))/dlog(qs(1))-1.e-16)/block_nx1)+1+nbefore
        !    if (x>=xprobmax1) igslice(level)=domain_nx1*2**(level-1)
        !  end do
        !else if(x .ge. -xmgrid(1,1)*half) then
        !  xnew=x+xmgrid(1,1)*half
        !  do level=1,refine_max_level
        !    nbefore=nstretchedblocks_baselevel(1)/2*2**(level-1)
        !    if (level .gt. 1) xlgrid(1,level)=xlgrid(1,level-1)/2.d0
 !igslice(level)=floor(((xnew-1.e-16)/xlgrid(1,level))/block_nx1)+1+nbefore
        !  end do
        !else
        !  xnew=xmgrid(1,1)*half-x
        !  do level=1,refine_max_level
        !    qs(1)=qstretch_baselevel(1)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(1,level)=xsgrid(1,level-1)/(one+qs(1))
 !igslice(level)=nstretchedblocks_baselevel(1)/2*2**(level-1) !      -floor((dlog(one-xnew/xsgrid(1,level)*(one-qs(1)))/dlog(qs(1)))/block_nx1)
        !    if (x<=xprobmin1) igslice(level)=1
        !  end do
        !end if

      case default
        call mpistop("stretch type not supported by get_igslice")
      end select

    case (2)
      select case (stretch_type(2))
      case (stretch_none) ! This dimension is unstretched
        if (x<=xprobmin2) then
          igslice=1
        else if (x>=xprobmax2) then
          do level=1,refine_max_level
            igslice(level)=ng2(level)
          end do
        else
          distance=x-xprobmin2
          do level = 1, refine_max_level
            igslice(level) = ceiling(distance/dg2(level))
          end do
        end if

      case (stretch_uni) ! Uniform stretching

        if (x<=xprobmin2) then
          igslice=1
        else if (x>=xprobmax2) then
          do level=1,refine_max_level
            igslice(level)=ng2(level)
          end do
        else
          distance=x-xprobmin2
          do level=1,refine_max_level
            igslice(level)= ceiling(dlog(distance/dxfirst(level,&
               2)*(qstretch(level,2)-1.d0)+1.d0)/dlog(qstretch(level,&
               2))/dble(block_nx2))
          end do
        end if

      case (stretch_symm) ! Symmetric stretching
        ! symmetric stretch about 0.5*(xprobmin+xprobmax)
        if (x<=xprobmin2) then
          igslice=1
        else if (x>=xprobmax2) then
          do level=1,refine_max_level
            igslice(level)=ng2(level)
          end do
        else if(x<xprobmin2+xstretch2) then
          distance=xprobmin2+xstretch2-x
          ! stretch to left from xprobmin+xstretch
          do level=1,refine_max_level
            igslice(level) = nstretchedblocks(level,&
               2)/2-int(dlog(distance*(qstretch(level,2)-one)/dxfirst(level,&
               2)+one)/dlog(qstretch(level,2))/dble(block_nx2))
          end do
        else if(x>xprobmax2-xstretch2) then
          distance=x-xprobmax2+xstretch2
          ! stretch to right from xprobmax-xstretch
          do level=1,refine_max_level
            igslice(level) = ceiling(dlog(distance*(qstretch(level,&
               2)-one)/dxfirst(level,2)+one)/dlog(qstretch(level,&
               2))/dble(block_nx2))+ng2(level)-nstretchedblocks(level,2)/2
          end do
        else
          ! possible non-stretched central part
          distance=x-xprobmin2-xstretch2
          do level=1,refine_max_level
            igslice(level)=nstretchedblocks(level,&
               2)/2+ceiling(distance/dg2(level))
          end do
        end if
 !xsgrid(2,1)=half*(xprobmax2-xprobmin2)!        /(half*domain_nx2-half*nstretchedblocks_baselevel(2)*block_nx2 !        +(one-qstretch_baselevel(2)**(half*nstretchedblocks_baselevel(2)*block_nx2)) !        /(one-qstretch_baselevel(2)))
 !xmgrid(2,1)=xsgrid(2,1)*(domain_nx2-nstretchedblocks_baselevel(2)*block_nx2)
        !xlgrid(2,1)=xsgrid(2,1)

        !if (x .ge. xmgrid(2,1)*half) then
        !  xnew=x-xmgrid(2,1)*half
        !  do level=1,refine_max_level
        !    qs(2)=qstretch_baselevel(2)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(2,level)=xsgrid(2,level-1)/(one+qs(2))
 !nbefore=(domain_nx2/block_nx2-nstretchedblocks_baselevel(2)/2)*2**(level-1)
 !igslice(level)=floor((dlog(one-xnew/xsgrid(2,level)*(one-qs(2)))/dlog(qs(2))-1.e-16)/block_nx2)+1+nbefore
        !    if (x>=xprobmax2) igslice(level)=domain_nx2*2**(level-1)
        !  end do
        !else if(x .ge. -xmgrid(2,1)*half) then
        !  xnew=x+xmgrid(2,1)*half
        !  do level=1,refine_max_level
        !    nbefore=nstretchedblocks_baselevel(2)/2*2**(level-1)
        !    if (level .gt. 1) xlgrid(2,level)=xlgrid(2,level-1)/2.d0
 !igslice(level)=floor(((xnew-1.e-16)/xlgrid(2,level))/block_nx2)+1+nbefore
        !  end do
        !else
        !  xnew=xmgrid(2,1)*half-x
        !  do level=1,refine_max_level
        !    qs(2)=qstretch_baselevel(2)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(2,level)=xsgrid(2,level-1)/(one+qs(2))
 !igslice(level)=nstretchedblocks_baselevel(2)/2*2**(level-1) !      -floor((dlog(one-xnew/xsgrid(2,level)*(one-qs(2)))/dlog(qs(2)))/block_nx2)
        !    if (x<=xprobmin2) igslice(level)=1
        !  end do
        !end if

      case default
        call mpistop("stretch type not supported by get_igslice")
      end select

    case (3)
      select case (stretch_type(3))
      case (stretch_none) ! This dimension is unstretched
        if (x<=xprobmin3) then
          igslice=1
        else if (x>=xprobmax3) then
          do level=1,refine_max_level
            igslice(level)=ng3(level)
          end do
        else
          distance=x-xprobmin3
          do level = 1, refine_max_level
            igslice(level) = ceiling(distance/dg3(level))
          end do
        end if

      case (stretch_uni) ! Uniform stretching

        if (x<=xprobmin3) then
          igslice=1
        else if (x>=xprobmax3) then
          do level=1,refine_max_level
            igslice(level)=ng3(level)
          end do
        else
          distance=x-xprobmin3
          do level=1,refine_max_level
            igslice(level)= ceiling(dlog(distance/dxfirst(level,&
               3)*(qstretch(level,3)-1.d0)+1.d0)/dlog(qstretch(level,&
               3))/dble(block_nx3))
          end do
        end if

      case (stretch_symm) ! Symmetric stretching
        ! symmetric stretch about 0.5*(xprobmin+xprobmax)
        if (x<=xprobmin3) then
          igslice=1
        else if (x>=xprobmax3) then
          do level=1,refine_max_level
            igslice(level)=ng3(level)
          end do
        else if(x<xprobmin3+xstretch3) then
          distance=xprobmin3+xstretch3-x
          ! stretch to left from xprobmin+xstretch
          do level=1,refine_max_level
            igslice(level) = nstretchedblocks(level,&
               3)/2-int(dlog(distance*(qstretch(level,3)-one)/dxfirst(level,&
               3)+one)/dlog(qstretch(level,3))/dble(block_nx3))
          end do
        else if(x>xprobmax3-xstretch3) then
          distance=x-xprobmax3+xstretch3
          ! stretch to right from xprobmax-xstretch
          do level=1,refine_max_level
            igslice(level) = ceiling(dlog(distance*(qstretch(level,&
               3)-one)/dxfirst(level,3)+one)/dlog(qstretch(level,&
               3))/dble(block_nx3))+ng3(level)-nstretchedblocks(level,3)/2
          end do
        else
          ! possible non-stretched central part
          distance=x-xprobmin3-xstretch3
          do level=1,refine_max_level
            igslice(level)=nstretchedblocks(level,&
               3)/2+ceiling(distance/dg3(level))
          end do
        end if
 !xsgrid(3,1)=half*(xprobmax3-xprobmin3)!        /(half*domain_nx3-half*nstretchedblocks_baselevel(3)*block_nx3 !        +(one-qstretch_baselevel(3)**(half*nstretchedblocks_baselevel(3)*block_nx3)) !        /(one-qstretch_baselevel(3)))
 !xmgrid(3,1)=xsgrid(3,1)*(domain_nx3-nstretchedblocks_baselevel(3)*block_nx3)
        !xlgrid(3,1)=xsgrid(3,1)

        !if (x .ge. xmgrid(3,1)*half) then
        !  xnew=x-xmgrid(3,1)*half
        !  do level=1,refine_max_level
        !    qs(3)=qstretch_baselevel(3)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(3,level)=xsgrid(3,level-1)/(one+qs(3))
 !nbefore=(domain_nx3/block_nx3-nstretchedblocks_baselevel(3)/2)*2**(level-1)
 !igslice(level)=floor((dlog(one-xnew/xsgrid(3,level)*(one-qs(3)))/dlog(qs(3))-1.e-16)/block_nx3)+1+nbefore
        !    if (x>=xprobmax3) igslice(level)=domain_nx3*2**(level-1)
        !  end do
        !else if(x .ge. -xmgrid(3,1)*half) then
        !  xnew=x+xmgrid(3,1)*half
        !  do level=1,refine_max_level
        !    nbefore=nstretchedblocks_baselevel(3)/2*2**(level-1)
        !    if (level .gt. 1) xlgrid(3,level)=xlgrid(3,level-1)/2.d0
 !igslice(level)=floor(((xnew-1.e-16)/xlgrid(3,level))/block_nx3)+1+nbefore
        !  end do
        !else
        !  xnew=xmgrid(3,1)*half-x
        !  do level=1,refine_max_level
        !    qs(3)=qstretch_baselevel(3)**(one/2.d0**(level-1))
        !    if (level .gt. 1) xsgrid(3,level)=xsgrid(3,level-1)/(one+qs(3))
 !igslice(level)=nstretchedblocks_baselevel(3)/2*2**(level-1) !      -floor((dlog(one-xnew/xsgrid(3,level)*(one-qs(3)))/dlog(qs(3)))/block_nx3)
        !    if (x<=xprobmin3) igslice(level)=1
        !  end do
        !end if

      case default
        call mpistop("stretch type not supported by get_igslice")
      end select

    case default
      call mpistop("slice direction not clear in get_igslice")
    end select

  end subroutine get_igslice