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=8) :: answer, lg^d
    integer(kind=4), intent(in) :: ig^d,ndim
    integer(kind=4) :: i
 
    ! 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
 
  !> Cost-weighted SFC partition. Reduces the per-rank block_cost array
  !> into a global Morton-indexed cost, EWMA-blends with the persistent
  !> costlist, and partitions the Morton order so each rank's cumulative
  !> cost is balanced. Falls back to equal-block-count if all costs are
  !> zero (cold start before the first measured step).
  subroutine get_morton_range_costed
    use mod_forest
    use mod_global_parameters
 
    integer :: ipe, Morton_no, igrid, ix
    double precision :: cost_total, cost_target, cost_cum
    double precision :: maxcount, mincount, maxcost, meancost
    double precision, allocatable :: cost_local(:)
    integer :: nblocks_per(0:npe-1)
 
    if (allocated(sfc_to_igrid)) deallocate(sfc_to_igrid)
    {#IFDEF EVOLVINGBOUNDARY
    if (allocated(sfc_phybound)) deallocate(sfc_phybound)
    }
 
    ! 1. Each rank fills its measured per-step block cost into the
    !    Morton-indexed cost_local at its own slots, zeros elsewhere.
    !    Morton numbering is invariant across load_balance migration, so
    !    costlist values built up via EWMA below stay meaningful when
    !    blocks change rank. Refinement events insert/remove Morton
    !    indices; the EWMA recovers within ~5 cycles for affected entries.
    allocate(cost_local(nleafs))
    cost_local = 0.0d0
    do morton_no = 1, nleafs
      if (sfc(2,morton_no) == mype) then
        igrid = sfc(1,morton_no)
        if (igrid > 0 .and. igrid <= max_blocks) then
          cost_local(morton_no) = block_cost(igrid)
        end if
      end if
    end do
 
    ! 2. Reduce so every rank has the global per-Morton cost vector. Since
    !    each Morton index is owned by exactly one rank in the pre-balance
    !    state, the SUM is just the value contributed by that one rank.
    if (npe > 1) then
      call mpi_allreduce(mpi_in_place, cost_local, nleafs, &
                         mpi_double_precision, mpi_sum, icomm, ierrmpi)
    end if
 
    ! 3. EWMA blend the measurement into the persistent global costlist.
    !    Skip slots with zero measurement (block didn't run advect this
    !    cycle for some reason — keep the prior estimate). This blend at
    !    the GLOBAL Morton level (not per-igrid) is what makes the
    !    partition migration-safe.
    do morton_no = 1, nleafs
      if (cost_local(morton_no) > 0.0d0) then
        costlist(morton_no) = lb_alpha * costlist(morton_no) &
                            + (1.0d0 - lb_alpha) * cost_local(morton_no)
      end if
    end do
    deallocate(cost_local)
 
    cost_total = sum(costlist(1:nleafs))
    if (cost_total <= 0.0d0) then
      ! Fully cold: fall back to equal-block partition.
      call get_morton_range
      return
    end if
 
    ! 4. Greedy cumulative-cost cut along the Morton-sorted leaves. Aim for
    !    each rank to receive cost_total / npe of cumulative work; cut at
    !    the leaf where the running sum first exceeds the per-rank target.
    nblocks_per = 0
    ipe = 0
    cost_cum = 0.0d0
    morton_start(0) = 1
    do morton_no = 1, nleafs
      cost_cum = cost_cum + costlist(morton_no)
      cost_target = (dble(ipe) + 1.0d0) * cost_total / dble(npe)
      nblocks_per(ipe) = nblocks_per(ipe) + 1
      if (cost_cum >= cost_target .and. ipe < npe-1) then
        morton_stop(ipe) = morton_no
        ipe = ipe + 1
        morton_start(ipe) = morton_no + 1
      end if
    end do
    morton_stop(npe-1) = nleafs
 
    ! 5. Diagnostics (also drive the existing log columns Xload/Xmemory).
    maxcount = dble(maxval(nblocks_per))
    mincount = dble(max(1,minval(nblocks_per)))
    xmemory = maxcount / mincount
    maxcost  = 0.0d0
    do ix = 0, npe-1
      maxcost = max(maxcost, sum(costlist(morton_start(ix):morton_stop(ix))))
    end do
    meancost = cost_total / dble(npe)
    xload    = maxcost / max(meancost, tiny(1.0d0))
 
    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_costed
 
end module mod_space_filling_curve