tree_amg_aggregate.f90

Go to the documentation of this file.

! Copyright (c) 2024, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module tree_amg_aggregate
  use tree_amg
  use utils
  use num_types
  use comm
  use mesh, only : mesh_t
  use logger, only : neko_log, log_size
  implicit none
 
  type, public :: tamg_agg_monitor_t
    integer :: level
    integer :: num_dofs
    integer :: num_aggs
    integer :: phase1_naggs
    integer :: phase1_ndof_aggregated
    integer :: phase2_naggs
    integer :: phase2_ndof_aggregated
  end type tamg_agg_monitor_t
 
contains
 
  subroutine aggregate_finest_level(tamg, lx, ly, lz, ne)
    type(tamg_hierarchy_t), intent(inout) :: tamg
    integer, intent(in) :: lx, ly, lz, ne
    integer :: i, j, k, l, nl, nt
    integer :: lid, gid_ptr
    integer :: lvl_id
    lvl_id = 1
    nl = lx*ly*lz
    nt = nl*ne
    call tamg_lvl_init(tamg%lvl(lvl_id), lvl_id, ne, nt)
    gid_ptr = 1
    do l = 1, ne
      tamg%lvl(lvl_id)%nodes_ptr(l) = gid_ptr
      tamg%lvl(lvl_id)%nodes_gid(l) = l
      call tamg_node_init( tamg%lvl(lvl_id)%nodes(l), l, nl)
      lid = 0
      do k = 1, lz
        do j = 1, ly
          do i = 1, lx
            lid = lid + 1
            tamg%lvl(lvl_id)%nodes(l)%dofs(lid) = linear_index(i,j,k,l,lx,ly,lz)
 
            tamg%lvl(lvl_id)%nodes_dofs(gid_ptr) = linear_index(i,j,k,l,lx,ly,lz)
            tamg%lvl(lvl_id)%nodes_gids(gid_ptr) = l
            gid_ptr = gid_ptr + 1
          end do
        end do
      end do
    end do
    tamg%lvl(lvl_id)%nodes_ptr(ne+1) = gid_ptr
 
    call aggregation_monitor_finest(lvl_id,nt,ne)
 
  end subroutine aggregate_finest_level
 
  subroutine agg_greedy_first_pass(naggs, max_aggs, n_elements, &
      facet_neigh, offset_el, n_facet, is_aggregated, aggregate_size)
    integer, intent(inout):: naggs
    integer, intent(in) :: max_aggs, n_elements
    integer, intent(in) :: facet_neigh(:,:)
    integer, intent(in) :: offset_el, n_facet
    integer, intent(inout) :: is_aggregated(:)
    integer, allocatable, intent(inout) :: aggregate_size(:)
    integer, allocatable :: as_tmp(:)
    real(kind=dp) :: random_value
    integer :: i, side, nhbr
    logical :: no_nhbr_agg
 
!-----!    do while (naggs .le. max_aggs)
!-----!      call random_number(random_value)
!-----!      i = floor(random_value * n_elements + 1)
    do i = 1, n_elements
      if (is_aggregated(i) .eq. -1) then
        no_nhbr_agg = .true.
        do side = 1, n_facet
          nhbr = facet_neigh(side, i) - offset_el
          if ((nhbr .gt. 0).and.(nhbr .le. n_elements)) then
            if (is_aggregated(nhbr) .ne. -1) then
              no_nhbr_agg = .false.
            end if
          end if
        end do! side
 
        if (no_nhbr_agg) then
          naggs = naggs + 1
          is_aggregated(i) = naggs
          if(size(aggregate_size).lt.naggs) then
            allocate(as_tmp(naggs + 20))
            as_tmp(1:size(aggregate_size)) = aggregate_size
            call move_alloc(as_tmp, aggregate_size)
          end if
          aggregate_size(naggs) = 1
          do side = 1, n_facet
            nhbr = facet_neigh(side, i) - offset_el
            if ((nhbr .gt. 0).and.(nhbr .le. n_elements)) then
              if (is_aggregated(nhbr) .eq. -1) then
                is_aggregated(nhbr) = naggs
                aggregate_size(naggs) = aggregate_size(naggs) + 1
              end if
            end if
          end do! side
        end if! no_nhbr_agg
      end if! is_aggregated(i)
    end do
  end subroutine agg_greedy_first_pass
 
  subroutine agg_greedy_second_pass(naggs, max_aggs, n_elements, &
      facet_neigh, offset_el, n_facet, is_aggregated, aggregate_size)
    integer, intent(inout):: naggs
    integer, intent(in) :: max_aggs, n_elements
    integer, intent(in) :: facet_neigh(:,:)
    integer, intent(in) :: offset_el, n_facet
    integer, intent(inout) :: is_aggregated(:)
    integer, intent(inout) :: aggregate_size(:)
    integer :: i, side, nhbr
    integer :: tnt_agg, tst_agg, tnt_size, tst_size
 
    do i = 1, n_elements
      if (is_aggregated(i) .eq. -1) then
        tnt_agg = -1
        tnt_size = 999!TODO: replace with large number
        tst_agg = -1
        tst_size = 999!TODO: replace with large number
        do side = 1, n_facet
          nhbr = facet_neigh(side, i) - offset_el
          if ((nhbr .gt. 0).and.(nhbr .le. n_elements)) then
            if (is_aggregated(nhbr) .ne. -1) then
              tst_agg = is_aggregated(nhbr)
              tst_size = aggregate_size(tst_agg)
              if (tst_size .lt. tnt_size) then
                tnt_size = tst_size
                tnt_agg = tst_agg
              end if
            end if
          end if
        end do
 
        if (tnt_agg .ne. -1) then
          is_aggregated(i) = tnt_agg
          aggregate_size(tnt_agg) = aggregate_size(tnt_agg) + 1
        else
          naggs = naggs + 1
          if (naggs .gt. size(aggregate_size)) then!TODO: another movealoc here? Error? the max_aggs needs to change though...
            call neko_error("Creating too many aggregates... something might be wrong... try increasing max_aggs")
          end if
          is_aggregated(i) = naggs
          aggregate_size(naggs) = 1
          do side = 1, n_facet
            nhbr = facet_neigh(side, i) - offset_el
            if ((nhbr .gt. 0).and.(nhbr .le. n_elements)) then
              if (is_aggregated(nhbr) .eq. -1) then
                aggregate_size(naggs) = aggregate_size(naggs) + 1
                is_aggregated(nhbr) = naggs
              end if
            end if
          end do
        end if
 
      end if
    end do
  end subroutine agg_greedy_second_pass
 
  subroutine agg_fill_nhbr_info(agg_nhbr, n_agg_nhbr, n_elements, &
      facet_neigh, offset_el, n_facet, is_aggregated, aggregate_size)
    integer, intent(inout) :: agg_nhbr(:,:)
    integer, intent(inout) :: n_agg_nhbr
    integer, intent(in) :: n_elements
    integer, intent(in) :: facet_neigh(:,:)
    integer, intent(in) :: offset_el, n_facet
    integer, intent(inout) :: is_aggregated(:)
    integer, intent(inout) :: aggregate_size(:)
    integer :: i, j, side, nhbr, tst_agg, tnt_agg
    logical :: agg_added
 
    n_agg_nhbr = 0
 
    do i = 1, n_elements!TODO: this is the lazy expensive way...
      tnt_agg = is_aggregated(i)
      do side = 1, n_facet
        nhbr = facet_neigh(side,i) - offset_el
        if ((nhbr .gt. 0).and.(nhbr .le. n_elements)) then
          tst_agg = is_aggregated(nhbr)
          if (tst_agg .le. 0) call neko_error("Unaggregated element detected. We do not want to handle that here...")
          if (tst_agg .ne. tnt_agg) then
            agg_added = .false.
            do j = 1, 50!TODO: this hard-coded value
              if ((agg_nhbr(j,tnt_agg) .eq. tst_agg)) then
                agg_added = .true.
              else if ((agg_nhbr(j,tnt_agg).eq.-1).and.(.not.agg_added)) then
                agg_nhbr(j,tnt_agg) = tst_agg
                agg_added = .true.
                n_agg_nhbr = max(n_agg_nhbr, j)
              end if
            end do! j
            if (.not.agg_added) call neko_error("Aggregates have too many neighbors... probably. Or some other error.")
          end if
        end if
      end do! side
    end do! i
  end subroutine agg_fill_nhbr_info
 
  subroutine aggregate_greedy(tamg, lvl_id, max_aggs, facet_neigh, agg_nhbr)
    type(tamg_hierarchy_t), intent(inout) :: tamg
    integer, intent(in) :: lvl_id
    integer, intent(in) :: max_aggs
    integer, intent(in) :: facet_neigh(:,:)
    integer, intent(inout), allocatable :: agg_nhbr(:,:)
    integer, allocatable :: is_aggregated(:)
    integer, allocatable :: aggregate_size(:)
    integer :: n_elements, naggs, n_facet, offset_el
    integer :: i, j, l, ntot, n_agg_facet, gid_ptr
 
    if (lvl_id .lt. 2) then
      call neko_error("For now, can only use greedy agg after elms have been aggregated to points (level 1)")
    else if (lvl_id .eq. 2) then
      n_facet = 6
      offset_el = tamg%msh%offset_el
    else
      n_facet = 50!TODO: this hard-coded value. how many neighbors can there be?
      offset_el = 0
    end if
 
    naggs = 0
    n_elements = tamg%lvl(lvl_id-1)%nnodes
    allocate( is_aggregated( n_elements ) )
    allocate( aggregate_size( max_aggs ) )
 
    is_aggregated = -1
 
    aggregate_size = huge(i)!999999
 
    call agg_greedy_first_pass(naggs, max_aggs, n_elements, &
      facet_neigh, offset_el, n_facet, &
      is_aggregated, aggregate_size)
 
    call aggregation_monitor_phase1(lvl_id, n_elements, naggs, is_aggregated)
 
    call agg_greedy_second_pass(naggs, max_aggs, n_elements, &
      facet_neigh, offset_el, n_facet, &
      is_aggregated, aggregate_size)
 
    call aggregation_monitor_phase2(lvl_id, n_elements, naggs, is_aggregated)
 
    if (.true.) then
      allocate( agg_nhbr(50, max_aggs*2) )!TODO: this hard-coded n_facet value (20)...
      agg_nhbr = -1
      call agg_fill_nhbr_info(agg_nhbr, n_agg_facet, n_elements, &
        facet_neigh, offset_el, n_facet, &
        is_aggregated, aggregate_size)
    end if
 
    ntot = 0
    do l = 1, naggs
      ntot = ntot + aggregate_size(l)
    end do
    call tamg_lvl_init( tamg%lvl(lvl_id), lvl_id, naggs, ntot)
    ntot = 0
    gid_ptr = 1
    do l = 1, naggs
      tamg%lvl(lvl_id)%nodes_ptr(l) = gid_ptr
      tamg%lvl(lvl_id)%nodes_gid(l) = l
      call tamg_node_init( tamg%lvl(lvl_id)%nodes(l), l, aggregate_size(l))
      j = 0
      do i = 1, n_elements!TODO: this is the lazy expensive way...
        if (is_aggregated(i) .eq. l) then
          j = j+1
          tamg%lvl(lvl_id)%nodes(l)%dofs(j) = i
 
          tamg%lvl(lvl_id)%nodes_dofs(gid_ptr) = i
          tamg%lvl(lvl_id)%nodes_gids(gid_ptr) = l
          gid_ptr = gid_ptr + 1
        end if
      end do
      if (j .ne. tamg%lvl(lvl_id)%nodes(l)%ndofs) then
        call neko_error("Aggregation problem. Not enough dofs in node.")
      end if
      ntot = ntot + aggregate_size(l)
    end do
    tamg%lvl(lvl_id)%nodes_ptr(naggs+1) = gid_ptr
 
    call aggregation_monitor_final(lvl_id,ntot,naggs)
 
    deallocate( is_aggregated )
    deallocate( aggregate_size )
  end subroutine aggregate_greedy
 
  subroutine aggregate_end( tamg, lvl_id)
    type(tamg_hierarchy_t), intent(inout) :: tamg
    integer, intent(in) :: lvl_id
    integer :: nt, i
    nt = tamg%lvl(lvl_id-1)%nnodes
    call tamg_lvl_init( tamg%lvl(lvl_id), lvl_id, 1, nt)
 
    call tamg_node_init( tamg%lvl(lvl_id)%nodes(1), 1, nt)
    do i = 1, tamg%lvl(lvl_id-1)%nnodes
      tamg%lvl(lvl_id)%nodes(1)%dofs(i) = tamg%lvl(lvl_id-1)%nodes(i)%gid
 
      tamg%lvl(lvl_id)%nodes_dofs(i) = tamg%lvl(lvl_id-1)%nodes(i)%gid
      tamg%lvl(lvl_id)%nodes_gids(i) = 1
    end do
 
    tamg%lvl(lvl_id)%nodes_ptr(1) = 1
    tamg%lvl(lvl_id)%nodes_ptr(2) = 2
    tamg%lvl(lvl_id)%nodes_gid(1) = 1
  end subroutine aggregate_end
 
  subroutine aggregation_monitor_finest(lvl,ndof,nagg)
    integer, intent(in) :: lvl,ndof,nagg
    integer :: na_max, na_min, na_avg, na_sum
    character(len=LOG_SIZE) :: log_buf
 
    write(log_buf, '(A8,I2,A37)') '-- level',lvl,'-- Aggregation: Element-as-Aggregate'
    !write(log_buf, '(A44)') 'Aggregation: Element-as-Aggregate'
    call neko_log%message(log_buf)
 
    call mpi_allreduce(nagg, na_max, 1, mpi_integer, mpi_max, neko_comm)
    call mpi_allreduce(nagg, na_min, 1, mpi_integer, mpi_min, neko_comm)
    call mpi_allreduce(nagg, na_sum, 1, mpi_integer, mpi_sum, neko_comm)
    na_avg = na_sum / pe_size
    write(log_buf, '(A35,I6,A1,I6,A1,I6,A1)') 'Number of Aggregates: (',na_min,',',na_avg,',',na_max,')'
    call neko_log%message(log_buf)
 
  end subroutine aggregation_monitor_finest
 
  subroutine aggregation_monitor_phase1(lvl,ndof,nagg,is_aggregated)
    integer, intent(in) :: lvl,ndof,nagg
    integer, intent(in) :: is_aggregated(:)
    integer :: num_aggregated, i, na_max, na_min, na_avg, na_sum
    real(kind=rp) :: agg_frac
    integer(kind=i8) :: glb_dof, glb_aggd, loc_dof, loc_aggd
    character(len=LOG_SIZE) :: log_buf
    num_aggregated = 0
    do i = 1, ndof
      if (is_aggregated(i) .gt. -1) then
        num_aggregated = num_aggregated + 1
      end if
    end do
 
    !write(log_buf, '(A8,I2,A24)') '-- level',lvl,'-- Aggregation: phase1'
    write(log_buf, '(A27)') 'Aggregation: phase1'
    call neko_log%message(log_buf)
 
    loc_aggd = int(num_aggregated, i8)
    loc_dof = int(ndof, i8)
    call mpi_allreduce(loc_dof, glb_dof, 1, mpi_integer8, mpi_sum, neko_comm)
    call mpi_allreduce(loc_aggd, glb_aggd, 1, mpi_integer8, mpi_sum, neko_comm)
    agg_frac = real(glb_aggd,rp) / real(glb_dof,rp)
 
    call mpi_allreduce(nagg, na_max, 1, mpi_integer, mpi_max, neko_comm)
    call mpi_allreduce(nagg, na_min, 1, mpi_integer, mpi_min, neko_comm)
    call mpi_allreduce(nagg, na_sum, 1, mpi_integer, mpi_sum, neko_comm)
    na_avg = na_sum / pe_size
    write(log_buf, '(A35,I6,A1,I6,A1,I6,A1)') 'Number of Aggregates: (',na_min,',',na_avg,',',na_max,')'
    call neko_log%message(log_buf)
 
    call mpi_allreduce(num_aggregated, na_max, 1, mpi_integer, mpi_max, neko_comm)
    call mpi_allreduce(num_aggregated, na_min, 1, mpi_integer, mpi_min, neko_comm)
    call mpi_allreduce(num_aggregated, na_sum, 1, mpi_integer, mpi_sum, neko_comm)
    na_avg = na_sum / pe_size
    write(log_buf, '(A35,I6,A1,I6,A1,I6,A1,F6.2)') 'Aggregated: (',na_min,',',na_avg,',',na_max,')',(agg_frac*100)
    call neko_log%message(log_buf)
  end subroutine aggregation_monitor_phase1
 
  subroutine aggregation_monitor_phase2(lvl,ndof,nagg,is_aggregated)
    integer, intent(in) :: lvl,ndof,nagg
    integer, intent(in) :: is_aggregated(:)
    integer :: num_aggregated, i, na_max, na_min, na_avg, na_sum
    real(kind=rp) :: agg_frac
    integer(kind=i8) :: glb_dof, glb_aggd, loc_dof, loc_aggd
    character(len=LOG_SIZE) :: log_buf
    num_aggregated = 0
    do i = 1, ndof
      if (is_aggregated(i) .gt. -1) then
        num_aggregated = num_aggregated + 1
      end if
    end do
    !write(log_buf, '(A8,I2,A24)') '-- level',lvl,'-- Aggregation: phase2'
    write(log_buf, '(A27)') 'Aggregation: phase2'
    call neko_log%message(log_buf)
 
    loc_aggd = int(num_aggregated, i8)
    loc_dof = int(ndof, i8)
    call mpi_allreduce(loc_dof, glb_dof, 1, mpi_integer8, mpi_sum, neko_comm)
    call mpi_allreduce(loc_aggd, glb_aggd, 1, mpi_integer8, mpi_sum, neko_comm)
    agg_frac = real(glb_aggd,rp) / real(glb_dof,rp)
 
    call mpi_allreduce(nagg, na_max, 1, mpi_integer, mpi_max, neko_comm)
    call mpi_allreduce(nagg, na_min, 1, mpi_integer, mpi_min, neko_comm)
    call mpi_allreduce(nagg, na_sum, 1, mpi_integer, mpi_sum, neko_comm)
    na_avg = na_sum / pe_size
    write(log_buf, '(A35,I6,A1,I6,A1,I6,A1)') 'Number of Aggregates: (',na_min,',',na_avg,',',na_max,')'
    call neko_log%message(log_buf)
 
    call mpi_allreduce(num_aggregated, na_max, 1, mpi_integer, mpi_max, neko_comm)
    call mpi_allreduce(num_aggregated, na_min, 1, mpi_integer, mpi_min, neko_comm)
    call mpi_allreduce(num_aggregated, na_sum, 1, mpi_integer, mpi_sum, neko_comm)
    na_avg = na_sum / pe_size
    write(log_buf, '(A35,I6,A1,I6,A1,I6,A1,F6.2)') 'Aggregated: (',na_min,',',na_avg,',',na_max,')',(agg_frac*100)
    call neko_log%message(log_buf)
  end subroutine aggregation_monitor_phase2
 
  subroutine aggregation_monitor_final(lvl,ndof,nagg)
    integer, intent(in) :: lvl,ndof,nagg
    character(len=LOG_SIZE) :: log_buf
    !TODO: calculate min and max agg size
    !write(log_buf, '(A8,I2,A23,I6)') '-- level',lvl,'-- Aggregation: Done.', nagg
    write(log_buf, '(A26,I6)') 'Aggregation: Done.', nagg
    call neko_log%message(log_buf)
  end subroutine aggregation_monitor_final
 
end module tree_amg_aggregate