mod_space_filling_curve.t

Go to the documentation of this file.

module mod_space_filling_curve
  implicit none
 
contains
 
  !> build Morton space filling curve for level 1 grid
  subroutine level1_morton_order
    ! use Morton curve to connect level 1 grid blocks
    use mod_forest
    use mod_global_parameters
 
    integer, allocatable :: gsq_sfc(:^D&),seq_sfc(:),seq_ig^D(:)
    integer :: ig^D, ngsq^D, isq, total_number
    !integer(kind=8), external :: mortonEncode
    logical, allocatable :: in_domain(:)
 
    ! use the smallest square/cube to cover the full domain 
    ngsq^d=2**ceiling(log(real(ng^d(1)))/log(2.0));
    {^nooned
    {ngsq^d=max(ngsq^dd) \}
    }
    total_number={ngsq^d|*}
    ! Morton number acquired by block numbers
    allocate(gsq_sfc(ngsq^d))
    ! Morton number in sequence
    allocate(seq_sfc(total_number))
    ! block numbers in sequence
    {allocate(seq_ig^d(total_number))\}
    allocate(in_domain(total_number))
    in_domain=.true.
    ! get Morton-order numbers in the square/cube
    {do ig^db=1,ngsq^db\}
       gsq_sfc(ig^d)=int(mortonencode(ig^d-1,ndim))+1
       seq_sfc(gsq_sfc(ig^d))=gsq_sfc(ig^d)
       {seq_ig^d(gsq_sfc(ig^dd))=ig^d \}
    {end do\}
    ! mark blocks that are out of the domain and change Morton number
    do isq=1,total_number
       if (seq_ig^d(isq)>ng^d(1)|.or.) then
         seq_sfc(isq:total_number)=seq_sfc(isq:total_number)-1
         in_domain(isq)=.false.
       end if
    end do
    ! copy the modified Morton numbers to the blocks in the domain
    if(.not. allocated(iglevel1_sfc)) allocate(iglevel1_sfc(ng^d(1)))
    if(.not. allocated(sfc_iglevel1)) allocate(sfc_iglevel1(ndim,nglev1))
    do isq=1,total_number
      if(in_domain(isq)) then
        iglevel1_sfc(seq_ig^d(isq))=seq_sfc(isq)
        {sfc_iglevel1(^d,seq_sfc(isq))=seq_ig^d(isq) \}
      end if
    end do
 
    deallocate(gsq_sfc,seq_sfc,seq_ig^d,in_domain)
 
  end subroutine level1_morton_order
 
  integer(kind=8) function mortonencode(ig^D,ndim)
    use iso_fortran_env, only : int64
    implicit none
    integer(kind=4), intent(in) :: ig^d,ndim
    integer(kind=4) :: i
    integer(kind=8) :: answer, lg^d
 
    ! Create a 64-bit version of ig^D
    lg^d=ig^d;
    answer = 0
 
    do i=0,64/ndim
      {^ifoned answer=ig1}
      {^iftwod
       answer=ior(answer,ior(ishft(iand(lg1,ishft(1_int64,i)),i),&
              ishft(iand(lg2,ishft(1_int64,i)),i+1)))
      \}
      {^ifthreed
       answer=ior(answer,ior(ishft(iand(lg1,ishft(1_int64,i)),2*i),ior(ishft(&
       iand(lg2,ishft(1_int64,i)),2*i+1),ishft(iand(lg3,ishft(1_int64,i)),2*i+2))))
      \}
    end do
    mortonencode=answer
    return
  end function mortonencode
 
  !> Construct Morton-order as a global recursive lexicographic ordering.
  subroutine amr_morton_order
    use mod_forest
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    integer :: ig^D, Morton_no, isfc
 
    morton_no=0
    nglev1={ng^d(1)*}
    do isfc=1,nglev1
       ig^d=sfc_iglevel1(^d,isfc)\ 
       call get_morton_number(tree_root(ig^d))
    end do
 
    if (morton_no/=nleafs) then
       call mpistop("error in amr_Morton_order: Morton_no/=nleafs")
    end if
 
    contains
 
    recursive subroutine get_morton_number(tree)
 
      type(tree_node_ptr) :: tree
 
      integer :: ic^d
 
      if (tree%node%leaf) then
         morton_no=morton_no+1
         sfc(1,morton_no)=tree%node%igrid
         sfc(2,morton_no)=tree%node%ipe
         if (tree%node%active) then 
            sfc(3,morton_no)=1 
         else 
            sfc(3,morton_no)=0 
         end if
         if(tree%node%ipe==mype) igrid_to_sfc(tree%node%igrid)=morton_no
      else
         {do ic^db=1,2\}
            call get_morton_number(tree%node%child(ic^d))
         {end do\}
      end if
 
    end subroutine get_morton_number
 
  end subroutine amr_morton_order
 
  !> Set the Morton range for each processor
  subroutine get_morton_range
    use mod_forest
    use mod_global_parameters
 
    integer :: ipe, blocks_left, procs_left, num_blocks
 
    if (allocated(sfc_to_igrid)) deallocate(sfc_to_igrid)
    {#IFDEF EVOLVINGBOUNDARY
    if (allocated(sfc_phybound)) deallocate(sfc_phybound)
    }
 
    blocks_left = nleafs
 
    do ipe = 0, npe-1
      if (ipe == 0) then
        morton_start(ipe) = 1
      else
        morton_start(ipe) = morton_stop(ipe-1) + 1
      end if
 
      ! Compute how many blocks this cpu should take
      procs_left = npe - ipe
      num_blocks = ceiling(blocks_left / dble(procs_left))
      morton_stop(ipe) = morton_start(ipe) + num_blocks - 1
      blocks_left = blocks_left - num_blocks
    end do
 
    allocate(sfc_to_igrid(morton_start(mype):morton_stop(mype)))
    {#IFDEF EVOLVINGBOUNDARY
    allocate(sfc_phybound(nleafs))
    sfc_phybound=0
    }
 
  end subroutine get_morton_range
 
  subroutine get_morton_range_active
    use mod_forest
    use mod_global_parameters
 
    ! Cut the sfc based on weighted decision.  
    ! Oliver Porth, 02.02.2012
 
    !!!Here you choose the weithts:!!
    integer, parameter :: wa=3, wp=1
    ! wp : Weight for passive block
    ! wa : Weight for active block
    ! wp=0 : balance load (active blocks) exactly and 
    ! don't care about memory imbalance
    ! wp=wa : balance memory exactly and don't care about load
    ! Maximum possible memory imbalance is X=wa/wp.
 
    ! If you run into memory issues, decrease this ratio.
    ! Scaling should be better if you allow for higher ratio.  
    ! Note also that passive cells still do regridding and boundary swap.  
    ! I have best results with a ratio 2:1, but it is problem dependent.
    ! It can't do magic though...
    ! Best to make sure that the sfc is properly aligned with your problem. 
 
    integer :: ipe, Morton_no
    integer :: Mtot, Mstop, Mcurr
    ! For debugging: 
    integer :: nactive(0:npe-1),npassive(0:npe-1)
    !double precision, save :: ptasum=0
 
    if (allocated(sfc_to_igrid)) deallocate(sfc_to_igrid)
    {#IFDEF EVOLVINGBOUNDARY
    if (allocated(sfc_phybound)) deallocate(sfc_phybound)
    }
 
    mtot  = nleafs_active*wa+(nleafs-nleafs_active)*wp
    ipe = 0 
    mcurr = 0
 
    nactive=0
    npassive=0
 
    morton_start(0) = 1
    do morton_no=1,nleafs
       ! This is where we ideally would like to make the cuts:
       mstop  = (ipe+1)*int(mtot/npe)+min(ipe+1,mod(mtot,npe))
       ! Build up mass:
       mcurr = mcurr + (wa*sfc(3,morton_no)+wp*(1-sfc(3,morton_no)))
 
       if (sfc(3,morton_no)==1) then 
          nactive(ipe) = nactive(ipe) +1
          else
             npassive(ipe) = npassive(ipe) +1
          end if
 
       if (mcurr >= mstop) then 
          morton_stop(ipe) = morton_no
          ipe = ipe +1
          if (ipe>=npe) exit
          morton_start(ipe) = morton_no + 1
       end if
    end do
 
    xmemory=dble(maxval(npassive+nactive))/&
         dble(minval(npassive+nactive))
    xload=dble(maxval(nactive))/&
         dble(minval(nactive))
 
    !ptasum = ptasum +dble(nleafs-nleafs_active)/dble(nleafs_active)
 
    !if (mype == 0) print*, 'nleafs_passive:',nleafs-nleafs_active, 'nleafs_active:',nleafs_active,'ratio:',dble(nleafs-nleafs_active)/dble(nleafs_active),'mean ratio:',ptasum/it
 
    if (morton_stop(mype)>=morton_start(mype)) then
       allocate(sfc_to_igrid(morton_start(mype):morton_stop(mype)))
    {#IFDEF EVOLVINGBOUNDARY
       allocate(sfc_phybound(nleafs))
       sfc_phybound=0
    }
    end if
 
  end subroutine get_morton_range_active
 
end module mod_space_filling_curve