buildClusterWalls.F90

Go to the documentation of this file.

subroutine buildclusterwalls(level, sps, useDual, walls, famList, nFamList)
 
    ! This routine will will build a global reduced surface mesh and ADT
    ! for each cluster. It can build using either the primal mesh or the
    ! dual mesh depending on the  useDual option.
 
    use adtbuild, only: buildserialquad
    use blockpointers
    use communication
    use inputphysics
    use inputtimespectral
    use oversetdata
    use inputoverset
    use utils, only: setpointers, echk, pointreduce
    use warping, only: getcgnsmeshindices
    use sorting, only: faminlist
    implicit none
 
    ! Input Variables
    integer(kind=intType), intent(in) :: level, sps, nFamList
    logical, intent(in) :: useDual
    type(oversetwall), intent(inout), dimension(nClusters), target :: walls
    integer(Kind=intType), intent(in) :: famList(nFamList)
 
    ! Local Variables
    integer(kind=intType) :: i, j, k, l, ii, jj, kk, nn, mm, iNode, iCell, c
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, ni, nj, nUnique, cellID, cellID2
    integer(kind=intType) :: ierr, iDim, lj
 
    ! Data for local surface
    integer(kind=intType) :: nNodes, nCells
    integer(kind=intType) :: nNodesLocal, nCellsLocal
    integer(kind=intType), dimension(:, :), allocatable :: connLocal
    integer(kind=intType), dimension(:), allocatable :: clusterNodeLocal
    integer(kind=intType), dimension(:), allocatable :: clusterCellLocal
    real(kind=realtype), dimension(:, :), allocatable :: nodeslocal
    real(kind=realtype), dimension(:, :, :), pointer :: xx, xx1, xx2, xx3, xx4 => null()
    integer(kind=intType), dimension(:, :, :), pointer :: globalCGNSNode => null()
    integer(kind=intType), dimension(:, :), pointer :: ind, indCGNS => null()
    integer(kind=intType), dimension(:, :), pointer :: indCell => null()
    logical :: regularOrdering
 
    ! Data for global surface
    integer(kind=intTYpe) :: nNodesGlobal, nCellsGlobal
    integer(kind=intType), dimension(:, :), allocatable, target :: connGlobal
    real(kind=realtype), dimension(:, :), allocatable, target :: nodesglobal
    integer(kind=intType), dimension(:), allocatable, target :: nodeIndicesGlobal
    integer(kind=intType), dimension(:), allocatable, target :: nodeIndicesCGNSGlobal
    integer(kind=intType), dimension(:), allocatable, target :: cellIndicesGlobal
 
    integer(kind=intType), dimension(:), allocatable :: nodesPerCluster, cellsPerCluster, cnc, ccc
    integer(kind=intType), dimension(:), allocatable :: clusterNodeGlobal
    integer(kind=intType), dimension(:), allocatable :: clusterCellGlobal
    integer(kind=intType), dimension(:), allocatable :: localNodeNums
    integer(kind=intType), dimension(:), allocatable :: nodeIndicesLocal
    integer(kind=intType), dimension(:), allocatable :: nodeIndicesCGNSLocal
    integer(kind=intType), dimension(:), allocatable :: cellIndicesLocal
    integer(kind=intType), dimension(:), allocatable :: cgnsIndices, curCGNSNode
 
    integer(kind=intType), dimension(:), allocatable :: nCellProc, cumCellProc
    integer(kind=intType), dimension(:), allocatable :: nNodeProc, cumNodeProc
    real(kind=realtype), dimension(:, :), allocatable :: uniquenodes
    integer(kind=intType), dimension(:), allocatable :: link
    real(kind=realtype), parameter :: tol = 1e-12
 
    ! Pointers for easier readibility
    integer(kind=intType), dimension(:, :), pointer :: conn
    integer(kind=intType), dimension(:), pointer :: tmpInd
 
    ! The first thing we do is gather all the surface nodes to
    ! each processor such that every processor can make it's own copy of
    ! the complete surface mesh to use to search. Note that this
    ! procedure *DOES NOT SCALE IN MEMORY*...ie eventually the surface
    ! mesh will become too large to store on a single processor,
    ! although this will probably not happen until the sizes get up in
    ! the hundreds of millions of cells.
 
    nnodeslocal = 0
    ncellslocal = 0
 
    ! Before we start generate a local node indices for the globalCGNS index
    ii = 0
    do nn = 1, ndom
        call setpointers(nn, level, sps)
        allocate (flowdoms(nn, level, sps)%globalCGNSNode(1:il, 1:jl, 1:kl))
        flowdoms(nn, level, sps)%globalCGNSNode = 0
        ii = ii + il * jl * kl
    end do
 
    if (level == 1) then
        allocate (cgnsindices(3 * ii))
        call getcgnsmeshindices(size(cgnsindices), cgnsindices)
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
                        ii = ii + 3
                        ! The reason for the +3 in the counter and the /3 is
                        ! that the CGNSMesh indicies include all DOF, so it
                        ! is the total number of mesh nodes *3. Here We only
                        ! care about nodes themselves so it is sufficient to
                        ! use basic ordering.
                        flowdoms(nn, level, sps)%globalCGNSNode(i, j, k) = cgnsindices(ii) / 3
                    end do
                end do
            end do
        end do
        deallocate (cgnsindices)
    end if
 
    do nn = 1, ndom
        call setpointers(nn, level, sps)
 
        do mm = 1, nbocos
            if (faminlist(bcdata(mm)%famID, famlist)) then
                ibeg = bcdata(mm)%inBeg
                iend = bcdata(mm)%inEnd
                jbeg = bcdata(mm)%jnBeg
                jend = bcdata(mm)%jnEnd
                if (usedual) then
                    nnodeslocal = nnodeslocal + &
                                  (iend - ibeg + 2) * (jend - jbeg + 2)
                    ncellslocal = ncellslocal + &
                                  (iend - ibeg + 1) * (jend - jbeg + 1)
                else
                    nnodeslocal = nnodeslocal + &
                                  (iend - ibeg + 1) * (jend - jbeg + 1)
                    ncellslocal = ncellslocal + &
                                  (iend - ibeg) * (jend - jbeg)
                end if
            end if
        end do
    end do
 
    ! Now communicate these sizes with everyone
    allocate (ncellproc(nproc), cumcellproc(0:nproc), &
              nnodeproc(nproc), cumnodeproc(0:nproc))
 
    call mpi_allgather(ncellslocal, 1, adflow_integer, ncellproc, 1, adflow_integer, &
                       adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgather(nnodeslocal, 1, adflow_integer, nnodeproc, 1, adflow_integer, &
                       adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Now make cumulative versions of these
    cumcellproc(0) = 0
    cumnodeproc(0) = 0
    do i = 1, nproc
        cumcellproc(i) = cumcellproc(i - 1) + ncellproc(i)
        cumnodeproc(i) = cumnodeproc(i - 1) + nnodeproc(i)
    end do
 
    ! And save the total number of nodes and cells for reference
    ncellsglobal = cumcellproc(nproc)
    nnodesglobal = cumnodeproc(nproc)
 
    ! Allocate the space for the local nodes and element connectivity
    allocate (nodeslocal(3, nnodeslocal), connlocal(4, ncellslocal), &
              clustercelllocal(ncellslocal), clusternodelocal(nnodeslocal), &
              nodeindiceslocal(nnodeslocal), nodeindicescgnslocal(nnodeslocal), &
              cellindiceslocal(ncellslocal))
 
    icell = 0
    inode = 0
    ! Second loop over the local walls
    do nn = 1, ndom
        call setpointers(nn, level, sps)
        c = clusters(cumdomproc(myid) + nn)
        globalcgnsnode => flowdoms(nn, level, sps)%globalCGNSNode
        do mm = 1, nbocos
            if (faminlist(bcdata(mm)%famID, famlist)) then
 
                jbeg = bcdata(mm)%jnBeg - 1; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg - 1; iend = bcdata(mm)%inEnd
 
                if (usedual) then
                    ! For the dual we have to allocate the pointer, xx.
                    select case (bcfaceid(mm))
                    case (imin)
                        xx1 => x(1, 0:jl, 0:kl, :)
                        xx2 => x(1, 1:je, 0:kl, :)
                        xx3 => x(1, 0:jl, 1:ke, :)
                        xx4 => x(1, 1:je, 1:ke, :)
                        ind => globalcell(2, 1:je, 1:ke)
                    case (imax)
                        xx1 => x(il, 0:jl, 0:kl, :)
                        xx2 => x(il, 1:je, 0:kl, :)
                        xx3 => x(il, 0:jl, 1:ke, :)
                        xx4 => x(il, 1:je, 1:ke, :)
                        ind => globalcell(il, 1:je, 1:ke)
 
                    case (jmin)
                        xx1 => x(0:il, 1, 0:kl, :)
                        xx2 => x(1:ie, 1, 0:kl, :)
                        xx3 => x(0:il, 1, 1:ke, :)
                        xx4 => x(1:ie, 1, 1:ke, :)
                        ind => globalcell(1:ie, 2, 1:ke)
 
                    case (jmax)
                        xx1 => x(0:il, jl, 0:kl, :)
                        xx2 => x(1:ie, jl, 0:kl, :)
                        xx3 => x(0:il, jl, 1:ke, :)
                        xx4 => x(1:ie, jl, 1:ke, :)
                        ind => globalcell(1:ie, jl, 1:ke)
 
                    case (kmin)
                        xx1 => x(0:il, 0:jl, 1, :)
                        xx2 => x(1:ie, 0:jl, 1, :)
                        xx3 => x(0:il, 1:je, 1, :)
                        xx4 => x(1:ie, 1:je, 1, :)
                        ind => globalcell(1:ie, 1:je, 2)
 
                    case (kmax)
                        xx1 => x(0:il, 0:jl, kl, :)
                        xx2 => x(1:ie, 0:jl, kl, :)
                        xx3 => x(0:il, 1:je, kl, :)
                        xx4 => x(1:ie, 1:je, kl, :)
                        ind => globalcell(1:ie, 1:je, kl)
                    end select
 
                else
                    select case (bcfaceid(mm))
                    case (imin)
                        xx => x(1, :, :, :)
                        ind => globalnode(1, :, :)
                        indcgns => globalcgnsnode(1, :, :)
                        ! Pointer to owned global cell indices
                        indcell => globalcell(2, :, :)
 
                    case (imax)
                        xx => x(il, :, :, :)
                        ind => globalnode(il, :, :)
                        indcgns => globalcgnsnode(il, :, :)
 
                        ! Pointer to owned global cell indices
                        indcell => globalcell(il, :, :)
 
                    case (jmin)
                        xx => x(:, 1, :, :)
                        ind => globalnode(:, 1, :)
                        indcgns => globalcgnsnode(:, 1, :)
                        ! Pointer to owned global cell indices
                        indcell => globalcell(:, 2, :)
 
                    case (jmax)
                        xx => x(:, jl, :, :)
                        ind => globalnode(:, jl, :)
                        indcgns => globalcgnsnode(:, jl, :)
                        ! Pointer to owned global cell indices
                        indcell => globalcell(:, jl, :)
 
                    case (kmin)
                        xx => x(:, :, 1, :)
                        ind => globalnode(:, :, 1)
                        indcgns => globalcgnsnode(:, :, 1)
                        ! Pointer to owned global cell indices
                        indcell => globalcell(:, :, 2)
 
                    case (kmax)
                        xx => x(:, :, kl, :)
                        ind => globalnode(:, :, kl)
                        indcgns => globalcgnsnode(:, :, kl)
 
                        ! Pointer to owned global cell indices
                        indcell => globalcell(:, :, kl)
 
                    end select
 
                    ! Just set hte 4 other pointers to xx so we can use the
                    ! same quarter-summation code below:
                    xx1 => xx
                    xx2 => xx
                    xx3 => xx
                    xx4 => xx
 
                end if
 
                ! We want to ensure that all the normals of the faces are
                ! consistent. To ensure this, we enforce that all normals
                ! are "into" the domain. Therefore we must treat difference
                ! faces of a block differently. For example for an iLow
                ! face, when looping over j-k in the regular way, results
                ! in in a domain inward pointing normal for iLow but
                ! outward pointing normal for iHigh. The same is true for
                ! kMin and kMax. However, it is reverse for the J-faces:
                ! This is becuase the way the pointers are extracted i then
                ! k is the reverse of what "should" be for consistency. The
                ! other two, the pointers are cyclic consistent: i,j->k,
                ! j,k (wrap) ->i, but for the j-direction is is i,k->j when
                ! to be consistent with the others it should be
                ! k,i->j. Hope that made sense.
 
                select case (bcfaceid(mm))
                case (imin, jmax, kmin)
                    regularordering = .true.
                case default
                    regularordering = .false.
                end select
 
                ! Now this can be reversed *again* if we have a block that
                ! is left handed.
                if (.not. righthanded) then
                    regularordering = .not. (regularordering)
                end if
 
                if (usedual) then
                    ! Start and end bounds for NODES
                    jbeg = bcdata(mm)%jnBeg - 1; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg - 1; iend = bcdata(mm)%inEnd
                else
                    ! Start and end bounds for NODES
                    jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                    ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
                end if
 
                ! ni, nj are the number of NODES
                ni = iend - ibeg + 1
                nj = jend - jbeg + 1
 
                ! Loop over the faces....this is the node sizes - 1
                if (regularordering) then
                    do j = 1, nj - 1
                        do i = 1, ni - 1
                            icell = icell + 1
                            connlocal(1, icell) = cumnodeproc(myid) + inode + (j - 1) * ni + i
                            connlocal(2, icell) = cumnodeproc(myid) + inode + (j - 1) * ni + i + 1
                            connlocal(3, icell) = cumnodeproc(myid) + inode + (j) * ni + i + 1
                            connlocal(4, icell) = cumnodeproc(myid) + inode + (j) * ni + i
                            ! Set the cluster
                            clustercelllocal(icell) = c
 
                            ! Save the global cell index
                            if (usedual) then
                                cellindiceslocal(icell) = 0
                            else
                                ! Valid only when using primary nodes
                                cellindiceslocal(icell) = indcell(ibeg + i + 1, jbeg + j + 1)
                            end if
                        end do
                    end do
                else
                    ! Do the reverse ordering
                    do j = 1, nj - 1
                        do i = 1, ni - 1
                            icell = icell + 1
                            connlocal(1, icell) = cumnodeproc(myid) + inode + (j - 1) * ni + i
                            connlocal(2, icell) = cumnodeproc(myid) + inode + (j) * ni + i
                            connlocal(3, icell) = cumnodeproc(myid) + inode + (j) * ni + i + 1
                            connlocal(4, icell) = cumnodeproc(myid) + inode + (j - 1) * ni + i + 1
 
                            ! Set the cluster
                            clustercelllocal(icell) = c
 
                            ! Save the global cell index
                            if (usedual) then
                                cellindiceslocal(icell) = 0
                            else
                                ! Valid only when using primary nodes
                                cellindiceslocal(icell) = indcell(ibeg + i + 1, jbeg + j + 1)
                            end if
                        end do
                    end do
                end if
 
                ! Loop over the nodes
                do j = jbeg, jend
                    do i = ibeg, iend
                        inode = inode + 1
                        ! The plus one is for the pointer offset
                        nodeslocal(:, inode) = fourth * ( &
                                               xx1(i + 1, j + 1, :) + xx2(i + 1, j + 1, :) + &
                                               xx3(i + 1, j + 1, :) + xx4(i + 1, j + 1, :))
 
                        clusternodelocal(inode) = c
                        nodeindiceslocal(inode) = ind(i + 1, j + 1) ! +1 for pointer offset
                        if (.not. usedual) then
                            nodeindicescgnslocal(inode) = indcgns(i, j) ! No pointer offset
                        else
                            nodeindicescgnslocal(inode) = 0
                        end if
                    end do
                end do
            end if
        end do
    end do
 
    ! Allocate space for the global reduced surface
    allocate (nodesglobal(3, nnodesglobal), connglobal(4, ncellsglobal), &
              clustercellglobal(ncellsglobal), clusternodeglobal(nnodesglobal), &
              nodeindicesglobal(nnodesglobal), nodeindicescgnsglobal(nnodesglobal), &
              cellindicesglobal(ncellsglobal))
 
    ! Communicate the nodes, connectivity and cluster information to everyone
    call mpi_allgatherv(nodeslocal, 3 * nnodeslocal, adflow_real, &
                        nodesglobal, nnodeproc * 3, cumnodeproc * 3, adflow_real, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgatherv(clusternodelocal, nnodeslocal, adflow_integer, &
                        clusternodeglobal, nnodeproc, cumnodeproc, adflow_integer, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgatherv(nodeindiceslocal, nnodeslocal, adflow_integer, &
                        nodeindicesglobal, nnodeproc, cumnodeproc, adflow_integer, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgatherv(nodeindicescgnslocal, nnodeslocal, adflow_integer, &
                        nodeindicescgnsglobal, nnodeproc, cumnodeproc, adflow_integer, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgatherv(connlocal, 4 * ncellslocal, adflow_integer, &
                        connglobal, ncellproc * 4, cumcellproc * 4, adflow_integer, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgatherv(clustercelllocal, ncellslocal, adflow_integer, &
                        clustercellglobal, ncellproc, cumcellproc, adflow_integer, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    call mpi_allgatherv(cellindiceslocal, ncellslocal, adflow_integer, &
                        cellindicesglobal, ncellproc, cumcellproc, adflow_integer, &
                        adflow_comm_world, ierr)
    call echk(ierr, __file__, __line__)
 
    ! Free the local data we do not need anymore
    deallocate (nodeslocal, connlocal, clustercelllocal, clusternodelocal, &
                ncellproc, cumcellproc, nnodeproc, cumnodeproc, nodeindiceslocal, &
                nodeindicescgnslocal, cellindiceslocal)
 
    ! We will now build separate trees for each cluster.
    allocate (nodespercluster(nclusters), cellspercluster(nclusters), &
              cnc(nclusters), ccc(nclusters))
    nodespercluster = 0
    cellspercluster = 0
 
    ! Count up the total number of elements and nodes on each cluster
    do i = 1, ncellsglobal
        cellspercluster(clustercellglobal(i)) = cellspercluster(clustercellglobal(i)) + 1
    end do
 
    do i = 1, nnodesglobal
        nodespercluster(clusternodeglobal(i)) = nodespercluster(clusternodeglobal(i)) + 1
    end do
 
    ! Create the list of the walls. We are reusing the overset wall derived type here.
    allocate (localnodenums(nnodesglobal))
 
    ! Allocate the memory for each of the cluster nodes
    do i = 1, nclusters
        nnodes = nodespercluster(i)
        ncells = cellspercluster(i)
        walls(i)%nCells = ncells
        walls(i)%nNodes = nnodes
 
        allocate (walls(i)%x(3, nnodes), walls(i)%conn(4, ncells), &
                  walls(i)%ind(nnodes))
        allocate (walls(i)%indCell(ncells))
    end do
 
    ! We now loop through the master list of nodes and elements and
    !  "push" them back to where they should go. We also keep track of
    !  the local node numbers so that the cluster surcells can update
    !  their own conn.
    localnodenums = 0
    cnc = 0
    do i = 1, nnodesglobal
        c = clusternodeglobal(i) ! Cluter this node belongs to
        cnc(c) = cnc(c) + 1 ! "cluster node count:" the 'nth' node for this  cluster
 
        walls(c)%x(:, cnc(c)) = nodesglobal(:, i)
        walls(c)%ind(cnc(c)) = nodeindicesglobal(i)
        localnodenums(i) = cnc(c)
    end do
 
    ccc = 0
    do i = 1, ncellsglobal
        c = clustercellglobal(i)
        ccc(c) = ccc(c) + 1 ! "Cluster cell count" the 'nth' cell for this cluster
        walls(c)%conn(:, ccc(c)) = connglobal(:, i)
 
        walls(c)%indCell(ccc(c)) = cellindicesglobal(i)
    end do
 
    do i = 1, nclusters
 
        ncells = walls(i)%nCells
        nnodes = walls(i)%nNodes
 
        ! Fistly we need to update the conn to use our local node ordering.
        do j = 1, ncells
            do k = 1, 4
                walls(i)%conn(k, j) = localnodenums(walls(i)%conn(k, j))
            end do
        end do
 
        ! Allocate temporary space for doing the point reduction.
        allocate (uniquenodes(3, nnodes), link(nnodes))
 
        call pointreduce(walls(i)%x, nnodes, tol, uniquenodes, link, nunique)
 
        ! Update the global indices. Use the returned link
        tmpind => walls(i)%ind
        allocate (walls(i)%ind(nunique))
        allocate (curcgnsnode(nunique))
        walls(i)%ind = -1
        curcgnsnode = -1
        do j = 1, walls(i)%nNodes
            ! Insted of blinding setting the index, we use the the
            ! nodeIndicesCGNSGlobal to only set the the globalNode with
            ! the smallest CGNS index. This guarantees that the same node
            ! ID is always selected independent of the block
            ! partitioning/splitting. Note that this will work even for
            ! the coarse levels when nodeIndicesCGNSGLobal are all 0's. In
            ! that case the first time wall(i)%ind(link(j)) is touched,
            ! that index is taken.
            lj = link(j)
            if (walls(i)%ind(lj) == -1) then
                ! Not set yet
                walls(i)%ind(lj) = tmpind(j)
                curcgnsnode(lj) = nodeindicescgnsglobal(j)
            else
                if (nodeindicescgnsglobal(j) < curcgnsnode(lj)) then
                    ! OR then potential global CGNS node index is LOWER
                    ! than the one I already have
                    walls(i)%ind(lj) = tmpind(j)
                    curcgnsnode(lj) = nodeindicescgnsglobal(j)
                end if
            end if
        end do
        deallocate (tmpind, curcgnsnode)
 
        ! Reset the number of nodes to be number of unique nodes
        nnodes = nunique
        walls(i)%nNodes = nnodes
 
        ! Update the nodes with the unique ones.
        do j = 1, nunique
            walls(i)%x(:, j) = uniquenodes(:, j)
        end do
 
        ! Update conn using the link:
        do j = 1, ncells
            do k = 1, 4
                walls(i)%conn(k, j) = link(walls(i)%conn(k, j))
            end do
        end do
 
        ! Unique nodes and link are no longer needed
        deallocate (link, uniquenodes)
 
        call buildserialquad(ncells, nnodes, walls(i)%x, walls(i)%conn, walls(i)%ADT)
    end do
 
    ! Clean up memeory
    deallocate (nodesglobal, connglobal, clustercellglobal, &
                clusternodeglobal, localnodenums, nodeindicesglobal, &
                nodeindicescgnsglobal)
 
    do nn = 1, ndom
        deallocate (flowdoms(nn, level, sps)%globalCGNSNode)
    end do
 
end subroutine buildclusterwalls