actuatorRegion.F90

Go to the documentation of this file.

module actuatorregion
 
    use constants
    use communication, only: commtype, internalcommtype
    use actuatorregiondata
    implicit none
 
contains
 
    subroutine computeactuatorregionvolume(nn, iRegion)
        use blockpointers, only: ndom, vol
        implicit none
 
        ! Inputs
        integer(kind=intType), intent(in) :: nn, iRegion
 
        ! Working
        integer(kind=intType) :: iii
        integer(kind=intType) :: i, j, k
 
        ! Loop over the region for this block
        do iii = actuatorregions(iregion)%blkPtr(nn - 1) + 1, actuatorregions(iregion)%blkPtr(nn)
            i = actuatorregions(iregion)%cellIDs(1, iii)
            j = actuatorregions(iregion)%cellIDs(2, iii)
            k = actuatorregions(iregion)%cellIDs(3, iii)
 
            ! Sum the volume of each cell within the region on this proc
            actuatorregions(iregion)%volLocal = actuatorregions(iregion)%volLocal + vol(i, j, k)
        end do
 
    end subroutine computeactuatorregionvolume
 
    ! ----------------------------------------------------------------------
    !                                                                      |
    !                    No Tapenade Routine below this line               |
    !                                                                      |
    ! ----------------------------------------------------------------------
 
#ifndef USE_TAPENADE
    subroutine addactuatorregion(pts, conn, axis1, axis2, famName, famID, &
                                 thrust, torque, heat, relaxStart, relaxEnd, nPts, nConn)
        ! Add a user-supplied integration surface.
 
        use communication, only: myid, adflow_comm_world
        use constants
        use adtbuild, only: buildserialquad, destroyserialquad
        use adtlocalsearch, only: mindistancetreesearchsinglepoint
        use adtutils, only: stack
        use adtdata
        use blockpointers, only: x, il, jl, kl, ndom, iblank, vol
        use adjointvars, only: ncellslocal
        use utils, only: setpointers, echk
        implicit none
 
        ! Input variables
        integer(kind=intType), intent(in) :: nPts, nConn, famID
        real(kind=realtype), dimension(3, nPts), intent(in), target :: pts
        integer(kind=intType), dimension(4, nConn), intent(in), target :: conn
        real(kind=realtype), intent(in), dimension(3) :: axis1, axis2
        character(len=*) :: famName
        real(kind=realtype) :: thrust, torque, heat, relaxstart, relaxend
 
        ! Working variables
        integer(kind=intType) :: i, j, k, nn, iDim, cellID, intInfo(3), sps, level, iii, ierr
        real(kind=realtype) :: dstar, frac, vollocal
        type(actuatorregiontype), pointer :: region
        real(kind=realtype), dimension(3) :: minx, maxx, v1, v2, v3, xcen, axisvec
        type(adttype) :: adt
        real(kind=realtype) :: axisvecnorm
        real(kind=realtype), dimension(:, :), allocatable :: norm
        integer(kind=intType), dimension(:), allocatable :: normcount
        integer(kind=intType), dimension(:, :), pointer :: tmp
 
        ! ADT Type required data
        integer(kind=intType), dimension(:), pointer :: frontleaves, frontleavesnew
        type(adtbboxtargettype), dimension(:), pointer :: BB
        real(kind=realtype) :: coor(4), uvw(5)
        real(kind=realtype) :: dummy(3, 2)
 
        nactuatorregions = nactuatorregions + 1
        if (nactuatorregions > nactuatorregionsmax) then
            print *, "Error: Exceeded the maximum number of actuatorDiskRegions.  "&
                 &"Increase nActuatorDiskRegionsMax"
            stop
        end if
 
        ! Save the input information
        region => actuatorregions(nactuatorregions)
        region%famName = famname
        region%famID = famid
        region%torque = torque
        region%heat = heat
        region%relaxStart = relaxstart
        region%relaxEnd = relaxend
        ! We use the axis to define the direction of F. Since we are
        ! dealing with rotating machinary, it is pretty good approximation
        ! to assume that the thrust is going to be in the direction of the
        ! axis.
        axisvec = axis2 - axis1
        axisvecnorm = sqrt((axisvec(1)**2 + axisvec(2)**2 + axisvec(3)**2))
        if (axisvecnorm < 1e-12) then
            print *, "Error: Axis cannot be determined by the supplied points. They are too close"
            stop
        end if
 
        axisvec = axisvec / axisvecnorm
 
        region%force = axisvec * thrust
        region%axisVec = axisvec
 
        allocate (region%blkPtr(0:ndom))
        region%blkPtr(0) = 0
 
        ! Next thing we need to do is to figure out if any of our cells
        ! are inside the actuator disk region. If so we will save them in
        ! the actuatorRegionType data structure
 
        ! Since this is effectively a wall-distance calc it gets super
        ! costly for the points far away. Luckly, we can do a fairly
        ! simple shortcut: Just compute the bounding box of the region and
        ! use that as the "already found" distance in the cloest point
        ! search. This will eliminate all the points further away
        ! immediately and this should be sufficiently fast.
 
        ! So...compute that bounding box:
        do idim = 1, 3
            minx(idim) = minval(pts(idim, :))
            maxx(idim) = maxval(pts(idim, :))
        end do
 
        ! Get the max distance. This should be quite conservative.
        dstar = (maxx(1) - minx(1))**2 + (maxx(2) - minx(2))**2 + (maxx(3) - minx(3))**2
 
        ! Now build the tree.
        call buildserialquad(size(conn, 2), size(pts, 2), pts, conn, adt)
 
        ! Compute the (averaged) unique nodal vectors:
        allocate (norm(3, size(pts, 2)), normcount(size(pts, 2)))
 
        norm = zero
        normcount = 0
 
        do i = 1, size(conn, 2)
 
            ! Compute cross product normal and normalize
            v1 = pts(:, conn(3, i)) - pts(:, conn(1, i))
            v2 = pts(:, conn(4, i)) - pts(:, conn(2, i))
 
            v3(1) = (v1(2) * v2(3) - v1(3) * v2(2))
            v3(2) = (v1(3) * v2(1) - v1(1) * v2(3))
            v3(3) = (v1(1) * v2(2) - v1(2) * v2(1))
            v3 = v3 / sqrt(v3(1)**2 + v3(2)**2 + v3(3)**2)
 
            ! Add to each of the four pts and increment the number added
            do j = 1, 4
                norm(:, conn(j, i)) = norm(:, conn(j, i)) + v3
                normcount(conn(j, i)) = normcount(conn(j, i)) + 1
            end do
        end do
 
        ! Now just divide by the norm count
        do i = 1, size(pts, 2)
            norm(:, i) = norm(:, i) / normcount(i)
        end do
 
        ! Norm count is no longer needed
        deallocate (normcount)
 
        ! Allocate the extra data the tree search requires.
        allocate (stack(100), bb(20), frontleaves(25), frontleavesnew(25))
 
        ! Allocate sufficient space for the maximum possible number of cellIDs
        allocate (region%cellIDs(3, ncellslocal(1)))
 
        ! Now search for all the coordinate. Note that We have explictly
        ! set sps to 1 becuase it is only implemented for single grid.
        sps = 1
        level = 1
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
            do k = 2, kl
                do j = 2, jl
                    do i = 2, il
                        ! Only check real cells
                        if (iblank(i, j, k) == 1) then
                            ! Compute the cell center
                            xcen = eighth * (x(i - 1, j - 1, k - 1, :) + x(i, j - 1, k - 1, :) &
                                             + x(i - 1, j, k - 1, :) + x(i, j, k - 1, :) &
                                             + x(i - 1, j - 1, k, :) + x(i, j - 1, k, :) &
                                             + x(i - 1, j, k, :) + x(i, j, k, :))
 
                            ! The current point to search for and continually
                            ! reset the "closest point already found" variable.
                            coor(1:3) = xcen
                            coor(4) = dstar
                            intinfo(3) = 0
                            call mindistancetreesearchsinglepoint(adt, coor, intinfo, &
                                                                  uvw, dummy, 0, bb, frontleaves, frontleavesnew)
                            cellid = intinfo(3)
                            if (cellid > 0) then
                                ! Now check if this was successful or now:
                                if (checkinside()) then
                                    ! Whoohoo! We are inside the region. Add this cell
                                    ! to the list.
                                    region%nCellIDs = region%nCellIDs + 1
                                    region%cellIDs(:, region%nCellIDs) = (/i, j, k/)
                                end if
                            end if
                        end if
                    end do
                end do
            end do
            ! Since we're doing all the blocks in order, simply store the
            ! current counter into blkPtr which gives up the range of cells
            ! we have found on this block
            region%blkPtr(nn) = region%nCellIDs
 
        end do
 
        ! Resize the cellIDs to the correct size now that we know the
        ! correct exact number.
        tmp => region%cellIDs
        allocate (region%cellIDs(3, region%nCellIDs))
        region%cellIDs = tmp(:, 1:region%nCellIDs)
        deallocate (tmp)
 
        ! Now go back and generate the total volume of the the cells we've flagged
        vollocal = zero
 
        do nn = 1, ndom
            call setpointers(nn, level, sps)
 
            ! Loop over the region for this block
            do iii = region%blkPtr(nn - 1) + 1, region%blkPtr(nn)
                i = region%cellIDs(1, iii)
                j = region%cellIDs(2, iii)
                k = region%cellIDs(3, iii)
                vollocal = vollocal + vol(i, j, k)
            end do
        end do
 
        call mpi_allreduce(vollocal, region%volume, 1, adflow_real, &
                           mpi_sum, adflow_comm_world, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Final memory cleanup
        deallocate (stack, norm, frontleaves, frontleavesnew, bb)
        call destroyserialquad(adt)
    contains
 
        function checkinside()
 
            implicit none
            logical :: checkinside
            integer(kind=intType) :: jj
            real(kind=realtype) :: shp(4), xp(3), normal(3), v1(3), dp
 
            ! bi-linear shape functions (CCW ordering)
            shp(1) = (one - uvw(1)) * (one - uvw(2))
            shp(2) = (uvw(1)) * (one - uvw(2))
            shp(3) = (uvw(1)) * (uvw(2))
            shp(4) = (one - uvw(1)) * (uvw(2))
 
            xp = zero
            normal = zero
            do jj = 1, 4
                xp = xp + shp(jj) * pts(:, conn(jj, cellid))
                normal = normal + shp(jj) * norm(:, conn(jj, cellid))
            end do
 
            ! Compute the dot product of normal with cell center
            ! (stored in coor) with the point on the surface.
            v1 = coor(1:3) - xp
            dp = normal(1) * v1(1) + normal(2) * v1(2) + normal(3) * v1(3)
 
            if (dp < zero) then
                checkinside = .true.
            else
                checkinside = .false.
            end if
        end function checkinside
    end subroutine addactuatorregion
 
    subroutine writeactuatorregions(fileName)
 
        ! This a (mostly) debug routine that is used to verify to the user
        ! the that the cells that the user thinks should be specified as
        ! being inside the actuator region actually are. We will dump a
        ! hex unstructured ascii tecplot file with all the zones we
        ! found. We won't be super concerned about efficiency here.
 
        use constants
        use utils, only: echk, pointreduce, setpointers
        use communication, only: myid, adflow_comm_world, nproc
        use blockpointers, only: x, ndom
        use commonformats, only: sci12
        implicit none
 
        ! Input
        character(len=*) :: fileName
 
        ! Working
        integer(kind=intType) :: iRegion, nn, i, j, k, ii, jj, kk, iii, kkk, iDim
        integer(kind=intType) :: level, sps, iProc, ierr, totalCount, offset, nUnique
        integer(kind=intType), dimension(:), allocatable :: sizesProc, cumSizesProc
        real(kind=realtype), dimension(:), allocatable :: pts, allpts
        real(kind=realtype), dimension(:, :), allocatable :: tmp, uniquepts
        real(kind=realtype), parameter :: tol = 1e-8
        integer(kind=intType), dimension(:), allocatable :: conn, allConn, link
        character(80) :: zoneName
        type(actuatorregiontype), pointer :: region
 
        ! Before we start the main region loop the root procesoor has to
        ! open up the tecplot file and write the header
 
        if (myid == 0) then
            open (unit=101, file=trim(filename), form='formatted')
            write (101, *) 'TITLE = "Actuator Regions"'
            write (101, *) 'Variables = "CoordinateX", "CoordinateY", "CoordinateZ"'
        end if
 
        ! Region Loop
        regionloop: do iregion = 1, nactuatorregions
 
            ! Only for the finest grid level.
            level = 1
            sps = 1
 
            ! Do an allgather with the number of actuator cells on each
            ! processor so that everyone knows the sizes and can compute the offsets,
            region => actuatorregions(iregion)
 
            allocate (sizesproc(nproc), cumsizesproc(0:nproc))
 
            call mpi_allgather(region%nCellIDs, 1, adflow_integer, sizesproc, 1, &
                               adflow_integer, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            cumsizesproc(0) = 0
            do iproc = 1, nproc
                cumsizesproc(iproc) = cumsizesproc(iproc - 1) + sizesproc(iproc)
            end do
 
            ! Fill up our own nodes/conn with the nodes we have here.
            allocate (conn(8 * region%nCellIDs), pts(24 * region%nCellIDs))
 
            kkk = 0
            do nn = 1, ndom
                call setpointers(nn, level, sps)
                ! Loop over the ranges for this block
                do iii = region%blkPtr(nn - 1) + 1, region%blkPtr(nn)
 
                    ! Carful with the conn values! They need to be in counter clock wise ordering!
                    offset = (iii - 1) * 8 + cumsizesproc(myid) * 8
                    conn((iii - 1) * 8 + 1) = 1 + offset
                    conn((iii - 1) * 8 + 2) = 2 + offset
                    conn((iii - 1) * 8 + 3) = 4 + offset
                    conn((iii - 1) * 8 + 4) = 3 + offset
                    conn((iii - 1) * 8 + 5) = 5 + offset
                    conn((iii - 1) * 8 + 6) = 6 + offset
                    conn((iii - 1) * 8 + 7) = 8 + offset
                    conn((iii - 1) * 8 + 8) = 7 + offset
 
                    ! Add in the 24 values for the nodal coordinates in coordinate
                    ! ordering. Do all the coordinates interlaced
                    do kk = -1, 0
                        do jj = -1, 0
                            do ii = -1, 0
                                do idim = 1, 3
                                    i = region%cellIDs(1, iii)
                                    j = region%cellIDs(2, iii)
                                    k = region%cellIDs(3, iii)
                                    kkk = kkk + 1
                                    pts(kkk) = x(i + ii, j + jj, k + kk, idim)
                                end do
                            end do
                        end do
                    end do
                end do
            end do
 
            ! Now that we've filled up our array, we can allocate the total
            ! space we need on the root proc and it
            if (myid == 0) then
                totalcount = sum(sizesproc)
                allocate (allconn(8 * totalcount), allpts(24 * totalcount))
            end if
 
            ! Perform the two gatherV's
            call mpi_gatherv(pts, region%nCellIDs * 24, adflow_real, &
                             allpts, 24 * sizesproc, 24 * cumsizesproc, adflow_real, &
                             0, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            call mpi_gatherv(conn, region%nCellIDs * 8, adflow_integer, &
                             allconn, 8 * sizesproc, 8 * cumsizesproc, adflow_integer, &
                             0, adflow_comm_world, ierr)
            call echk(ierr, __file__, __line__)
 
            ! We can deallocate all the per-proc memory now
            deallocate (sizesproc, cumsizesproc, pts, conn)
 
            if (myid == 0) then
 
                ! Now the poor root processor dumps everything out to a
                ! file. To help cut down on the already bloated file size,
                ! we'll point reduce it which will help tecplot display them
                ! better as well.
 
                allocate (tmp(3, totalcount * 8))
                do i = 1, totalcount * 8
                    do idim = 1, 3
                        tmp(idim, i) = allpts((i - 1) * 3 + idim)
                    end do
                end do
                deallocate (allpts)
                allocate (uniquepts(3, totalcount * 8), link(totalcount * 8))
 
                ! Get unique set of nodes.
                call pointreduce(tmp, totalcount * 8, tol, uniquepts, link, nunique)
 
                write (zonename, "(a,a,a)") 'Zone T="', trim(region%famName), ' Region"'
                write (101, *) trim(zonename)
                write (101, *) "Nodes = ", nunique, " Elements= ", totalcount, " ZONETYPE=FEBRICK"
                write (101, *) "DATAPACKING=BLOCK, VARLOCATION=([1,2,3]=NODAL, [4]=CELLCENTERED)"
 
                ! Write all the coordinates...this is horrendously slow...
                do idim = 1, 3
                    do i = 1, nunique
                        write (101, sci12) uniquepts(idim, i)
                    end do
                end do
 
                ! Write out the connectivity
                do i = 1, totalcount
                    do j = 1, 8
                        write (101, "(I8)", advance='no') link(allconn((i - 1) * 8 + j))
                    end do
                    write (101, "(1x)")
                end do
 
                ! Ditch the memory only allocated on this proc
                deallocate (allconn, link, tmp, uniquepts)
            end if
        end do regionloop
 
        ! Close the output file on the root proc
        if (myid == 0) then
            close (101)
        end if
    end subroutine writeactuatorregions
 
    subroutine integrateactuatorregions(localValues, famList, sps)
        !--------------------------------------------------------------
        ! Manual Differentiation Warning: Modifying this routine requires
        ! modifying the hand-written forward and reverse routines.
        ! --------------------------------------------------------------
 
        ! Perform volume integrals over the actuator region.
        use constants
        use blockpointers, only: vol, dw, w, ndom
        use flowvarrefstate, only: pref, uref
        use utils, only: setpointers
        use sorting, only: faminlist
        use residuals, only: sourceterms_block
        use actuatorregiondata
        use communication
        ! Input/output Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
 
        ! Working
        integer(kind=intType) :: nn, iRegion
        real(kind=realtype) :: plocal, plocald
        ! Zero the accumulation variable. We comptue flow power across
        ! 'all' actuaor zones. The famInclude is used to section out which
        ! one we want.
        plocal = zero
 
        domainloop: do nn = 1, ndom
            call setpointers(nn, 1, sps)
 
            ! Loop over each region
            regionloop: do iregion = 1, nactuatorregions
 
                ! Check if this needs to be included:
                faminclude: if (faminlist(actuatorregions(iregion)%famID, famlist)) then
 
                    ! If so, call the regular sourceTerms_block routine
                    call sourceterms_block(nn, .false., iregion, plocal)
 
                end if faminclude
            end do regionloop
        end do domainloop
 
        ! Add in the contribution from this processor.
        localvalues(ipower) = localvalues(ipower) + plocal
 
    end subroutine integrateactuatorregions
#ifndef USE_COMPLEX
    subroutine integrateactuatorregions_d(localValues, localValuesd, famList, sps)
        !--------------------------------------------------------------
        ! Manual Differentiation Warning: Modifying this routine requires
        ! modifying the hand-written forward and reverse routines.
        ! --------------------------------------------------------------
 
        ! Perform volume integrals over the actuator region.
        use constants
        use blockpointers, only: vol, dw, w, ndom
        use flowvarrefstate, only: pref, uref
        use utils, only: setpointers_d
        use sorting, only: faminlist
        use actuatorregiondata
        use residuals_d, only: sourceterms_block_d
 
        ! Input/output Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
 
        ! Working
        integer(kind=intType) :: nn, iRegion
        real(kind=realtype) :: plocal, plocald
 
        ! Zero the accumulation variable. We comptue flow power across
        ! 'all' actuaor zones. The famInclude is used to section out which
        ! one we want.
        plocal = zero
        plocald = zero
 
        domainloop: do nn = 1, ndom
            call setpointers_d(nn, 1, sps)
 
            ! Loop over each region
            regionloop: do iregion = 1, nactuatorregions
 
                ! Check if this needs to be included:
                faminclude: if (faminlist(actuatorregions(iregion)%famID, famlist)) then
 
                    ! If so, call the regular sourceTerms_block routine
                    call sourceterms_block_d(nn, .false., iregion, plocal, plocald)
 
                end if faminclude
            end do regionloop
        end do domainloop
 
        ! Add in the contribution from this processor.
        localvalues(ipower) = localvalues(ipower) + plocal
        localvaluesd(ipower) = localvaluesd(ipower) + plocald
 
    end subroutine integrateactuatorregions_d
 
    subroutine integrateactuatorregions_b(localValues, localValuesd, famList, sps)
        !--------------------------------------------------------------
        ! Manual Differentiation Warning: Modifying this routine requires
        ! modifying the hand-written forward and reverse routines.
        ! --------------------------------------------------------------
 
        ! Perform volume integrals over the actuator region.
        use constants
        use blockpointers, only: vol, dw, w, ndom
        use flowvarrefstate, only: pref, uref
        use utils, only: setpointers_b
        use sorting, only: faminlist
        use actuatorregiondata
        use residuals_b, only: sourceterms_block_b
 
        ! Input/output Variables
        real(kind=realtype), dimension(nLocalValues), intent(inout) :: localvalues, localvaluesd
        integer(kind=intType), dimension(:), intent(in) :: famList
        integer(kind=intType), intent(in) :: sps
 
        ! Working
        integer(kind=intType) :: nn, iRegion
        real(kind=realtype) :: plocal, plocald
 
        ! Zero the accumulation variable. We comptue flow power across
        ! 'all' actuaor zones. The famInclude is used to section out which
        ! one we want.
        plocal = zero
        plocald = zero
 
        ! Pull out the seed
        plocald = localvaluesd(ipower)
 
        domainloop: do nn = 1, ndom
            call setpointers_b(nn, 1, sps)
 
            ! Loop over each region
            regionloop: do iregion = 1, nactuatorregions
 
                ! Check if this needs to be included:
                faminclude: if (faminlist(actuatorregions(iregion)%famID, famlist)) then
 
                    ! If so, call the regular sourceTerms_block routine
                    call sourceterms_block_b(nn, .false., iregion, plocal, plocald)
 
                end if faminclude
            end do regionloop
        end do domainloop
    end subroutine integrateactuatorregions_b
 
#endif
#endif
end module actuatorregion