warping.F90

Go to the documentation of this file.

module warping
 
    ! This module cotains the required inferface functions for using an
    ! external mesh warping utility with ADflow
 
contains
 
    subroutine getcgnsmeshindices(ndof, indices)
 
        use constants
        use blockpointers, only: ndom, nbkglobal, il, jl, kl, ibegor, jbegor, kbegor
        use cgnsgrid, only: cgnsdoms, cgnsndom
        use utils, only: setpointers
        implicit none
 
        ! subroutine arguments
        integer(kind=intType), intent(in) :: ndof
        integer(kind=intType), intent(out) :: indices(ndof)
 
        ! Local Variables
        integer(kind=intType) :: nn, i, j, k, ii, indx, indy, indz, il_cg, jl_cg, kl_cg
        integer(kind=intType), allocatable, dimension(:) :: dof_offset
 
        allocate (dof_offset(cgnsndom))
        dof_offset(1) = 0
        do nn = 2, cgnsndom
            dof_offset(nn) = dof_offset(nn - 1) + &
                             cgnsdoms(nn - 1)%il * cgnsdoms(nn - 1)%jl * cgnsdoms(nn - 1)%kl * 3
        end do
 
        ii = 0
        do nn = 1, ndom
            call setpointers(nn, 1_inttype, 1_inttype)
 
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
                        il_cg = cgnsdoms(nbkglobal)%il
                        jl_cg = cgnsdoms(nbkglobal)%jl
                        kl_cg = cgnsdoms(nbkglobal)%kl
 
                        ii = ii + 1
 
                        indx = ibegor + i - 1
                        indy = jbegor + j - 1
                        indz = kbegor + k - 1
 
                        indices(ii * 3 - 2) = dof_offset(nbkglobal) + &
                                              (indz - 1) * jl_cg * il_cg * 3 + &
                                              (indy - 1) * il_cg * 3 + &
                                              (indx - 1) * 3
 
                        indices(ii * 3 - 1) = dof_offset(nbkglobal) + &
                                              (indz - 1) * jl_cg * il_cg * 3 + &
                                              (indy - 1) * il_cg * 3 + &
                                              (indx - 1) * 3 + 1
 
                        indices(ii * 3) = dof_offset(nbkglobal) + &
                                          (indz - 1) * jl_cg * il_cg * 3 + &
                                          (indy - 1) * il_cg * 3 + &
                                          (indx - 1) * 3 + 2
 
                    end do ! i loop
                end do ! j loop
            end do ! k loop
        end do ! domain loop
        deallocate (dof_offset)
    end subroutine getcgnsmeshindices
 
    subroutine setgrid(grid, ndof)
 
        ! The purpose of this routine is to set the grid dof as returned by
        ! the external warping. This function takes the "Base" grid at the
        ! first time instance and does rotation/translation operations to
        ! get the grid at subsequent time instances
        use constants
        use blockpointers, only: ndom, il, jl, kl, x
        use inputtimespectral, only: ntimeintervalsspectral
        use section, only: sections, nsections
        use inputtimespectral, only: ntimeintervalsspectral
        use monitor, only: timeunsteadyrestart, timeunsteady
        use inputphysics, only: equationmode
        use utils, only: setpointers, rotmatrixrigidbody
        use preprocessingapi, only: xhalo
        implicit none
 
        integer(kind=intType), intent(in) :: ndof
        real(kind=realtype), intent(in) :: grid(ndof)
 
        ! Local Variables
 
        integer(kind=intType) :: nn, i, j, k, counter, sps
        real(kind=realtype) :: t(nsections), dt(nsections)
        real(kind=realtype) :: displ(3)
        real(kind=realtype) :: told, tnew
 
        real(kind=realtype), dimension(3) :: rotationpoint, r
        real(kind=realtype), dimension(3, 3) :: rotationmatrix
 
        if (equationmode == steady .or. equationmode == timespectral) then
            timeunsteady = zero
 
            ! This is very straight forward...loop over all domains and set all elements
            do nn = 1, nsections
                dt(nn) = sections(nn)%timePeriod &
                         / real(ntimeintervalsspectral, realtype)
            end do
 
            do sps = 1, ntimeintervalsspectral
                do nn = 1, nsections
                    t(nn) = (sps - 1) * dt(nn)
                end do
 
                ! Compute the displacements due to the rigid motion of the mesh.
 
                displ(:) = zero
 
                tnew = timeunsteady + timeunsteadyrestart
                told = tnew - t(1)
 
                call rotmatrixrigidbody(tnew, told, rotationmatrix, rotationpoint)
                counter = 0
                do nn = 1, ndom
                    call setpointers(nn, 1_inttype, sps)
                    do k = 1, kl
                        do j = 1, jl
                            do i = 1, il
                                ! r is distance from grid point to rotationPoint
                                r = grid(3 * counter + 1:3 * counter + 3) - rotationpoint
 
                                x(i, j, k, :) = rotationpoint + matmul(rotationmatrix, r) + displ
                                counter = counter + 1
 
                            end do
                        end do
                    end do
                end do
                call xhalo(1_inttype)
            end do
        else
            counter = 0
            sps = 1
            do nn = 1, ndom
                call setpointers(nn, 1_inttype, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            x(i, j, k, :) = grid(3 * counter + 1:3 * counter + 3)
                            counter = counter + 1
                        end do
                    end do
                end do
            end do
            call xhalo(1_inttype)
        end if
 
    end subroutine setgrid
 
    subroutine setgridforoneinstance(grid, sps)
 
        ! The purpose of this routine is to set the grid dof as returned by
        ! the external warping. This routine will take in the deformed mesh
        ! and set it to "sps"th time instance
        use constants
        use blockpointers, only: ndom, il, jl, kl, x
        use section, only: sections, nsections
        use inputphysics, only: equationmode
        use utils, only: setpointers
        use preprocessingapi, only: xhalo
        implicit none
 
        real(kind=realtype), dimension(:), intent(in) :: grid
        integer, intent(in) :: sps
 
        ! Local Variables
 
        integer(kind=intType) :: nn, i, j, k, counter
 
        counter = 0
        do nn = 1, ndom
            call setpointers(nn, 1_inttype, sps)
            do k = 1, kl
                do j = 1, jl
                    do i = 1, il
                        x(i, j, k, :) = grid(3 * counter + 1:3 * counter + 3)
                        counter = counter + 1
                    end do
                end do
            end do
        end do
        call xhalo(1_inttype)
 
    end subroutine setgridforoneinstance
 
    subroutine getgrid(grid, ndof)
 
        ! Opposite of setGrid. This is ONLY a debugging function. NOT used
        ! in regular usage. Really only useful for direct mesh manipulation
        ! on single block and a single processor. s
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, x
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers
        implicit none
        integer(kind=intType), intent(in) :: ndof
        real(kind=realtype), intent(out) :: grid(ndof)
 
        ! Local Variables
        integer(kind=intType) :: nn, i, j, k, l, counter, sps
 
        ! This is very straight forward...loop over all domains and copy out
        counter = 1
        do sps = 1, ntimeintervalsspectral
            do nn = 1, ndom
                call setpointers(nn, 1_inttype, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                grid(counter) = x(i, j, k, l)
                                counter = counter + 1
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
    end subroutine getgrid
 
    subroutine getstateperturbation(randVec, nRand, randState, nRandState)
 
        use constants
        use cgnsgrid, only: cgnsdoms, cgnsndom
        use blockpointers, only: ndom, il, jl, kl, nx, ny, nz, x, nbkglobal, ibegor, jbegor, kbegor
        use communication, only: adflow_comm_world, myid
        use inputtimespectral, only: ntimeintervalsspectral
        use adjointvars, only: ncellslocal
        use flowvarrefstate, only: nw
        use utils, only: setpointers, echk
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: nRand, nRandState
        real(kind=realtype), intent(in), dimension(nRand) :: randvec
 
        ! Ouput Parameters
        real(kind=realtype), intent(out), dimension(nRandState) :: randstate
 
        ! Working parameters
        integer(kind=intType) :: i, j, k, ierr, l, nx_cg, ny_cg, nz_cg
        integer(kind=intType) :: sps, ii, indx, indy, indz, nn, cgnsInd
        integer(kind=intType) :: dofCGNSPerInstance
        integer(kind=intType), allocatable, dimension(:) :: dof_offset
 
        allocate (dof_offset(cgnsndom))
        dof_offset(1) = 0
        do nn = 2, cgnsndom
            dof_offset(nn) = dof_offset(nn - 1) + &
                             cgnsdoms(nn - 1)%nx * cgnsdoms(nn - 1)%ny * cgnsdoms(nn - 1)%nz * nw
        end do
 
        dofcgnsperinstance = nrand / ntimeintervalsspectral
 
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, 1, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
 
                            nx_cg = cgnsdoms(nbkglobal)%nx
                            ny_cg = cgnsdoms(nbkglobal)%ny
                            nz_cg = cgnsdoms(nbkglobal)%nz
 
                            indx = ibegor + i - 2
                            indy = jbegor + j - 2
                            indz = kbegor + k - 2
 
                            do l = 1, nw
                                cgnsind = (sps - 1) * dofcgnsperinstance + &
                                          dof_offset(nbkglobal) + &
                                          (indz - 1) * ny_cg * nx_cg * nw + &
                                          (indy - 1) * nx_cg * nw + &
                                          (indx - 1) * nw + l
                                randstate(nw * ii + l) = randvec(cgnsind)
                            end do
 
                            ii = ii + 1
 
                        end do
                    end do
                end do
            end do
        end do
    end subroutine getstateperturbation
 
    subroutine getsurfaceperturbation(xRand, nRand, randSurface, nRandSurface, famList, nFamList, sps)
 
        use constants
        use blockpointers, only: ndom, bcdata, nbocos, bcfaceid, il, jl, kl
        use communication, only: adflow_comm_world, myid
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers, echk
        use sorting, only: faminlist
        use oversetdata, only: zippermeshes, zippermesh, oversetpresent
        use surfacefamilies, only: bcfamgroups, familyexchange, bcfamexchange
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input Parameters
        integer(kind=intType), intent(in) :: nRand, nRandSurface
        real(kind=realtype), intent(in), dimension(nRand) :: xrand
        integer(kind=intType), intent(in) :: nFamList, famList(nFamList), sps
        ! Ouput Parameters
        real(kind=realtype), intent(inout), dimension(3*nRandSurface) :: randsurface
 
        ! Working parameters
        integer(kind=intType) :: i, j, k, ierr, iDim, iBeg, iEnd, jBeg, jEnd, nn, mm
        integer(kind=intType) :: ii, jj, indI, indJ, indK, jjInd, iBCGroup
        type(zippermesh), pointer :: zipper
        type(familyexchange), pointer :: exch
        logical :: BCGroupNeeded
        real(kind=realtype), dimension(:), pointer :: localptr
 
        ii = 0
        jj = 0
        domains: do nn = 1, ndom
            call setpointers(nn, 1_inttype, sps)
 
            ! Loop over the number of boundary subfaces of this block.
            bocos: do mm = 1, nbocos
 
                ! NODE Based
                jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
 
                faminclude: if (faminlist(bcdata(mm)%famID, famlist)) then
 
                    do j = jbeg, jend ! This is a node loop
                        do i = ibeg, iend ! This is a node loop
                            select case (bcfaceid(mm))
                            case (imin)
                                indi = 1
                                indj = i
                                indk = j
                            case (imax)
                                indi = il
                                indj = i
                                indk = j
                            case (jmin)
                                indi = i
                                indj = 1
                                indk = j
                            case (jmax)
                                indi = i
                                indj = jl
                                indk = j
                            case (kmin)
                                indi = i
                                indj = j
                                indk = 1
                            case (kmax)
                                indi = i
                                indj = j
                                indk = kl
                            end select
 
                            do idim = 1, 3
                                jjind = jj + (indk - 1) * il * jl + (indj - 1) * il + (indi - 1) + idim
                                randsurface(3 * ii + idim) = xrand(jjind)
                            end do
                            ii = ii + 1
                        end do
                    end do
                end if faminclude
            end do bocos
 
            ! jj is the counter through xRand. Increment it by the full
            ! block.
            jj = jj + il * jl * kl * 3
        end do domains
 
        ! No overset or not zipper, return
        if (.not. oversetpresent) then ! .or. .not. includeZipper) then
            return
        end if
 
        ! If there are zipper meshes, we must include the nodes that the
        ! zipper triangles will use.
        do ibcgroup = 1, nfamexchange
 
            zipper => zippermeshes(ibcgroup)
 
            if (.not. zipper%allocated) then
                cycle
            end if
 
            exch => bcfamexchange(ibcgroup, sps)
            bcgroupneeded = .false.
            bcgroupfamloop: do i = 1, size(bcfamgroups(ibcgroup)%famList)
                if (faminlist(bcfamgroups(ibcgroup)%famList(i), famlist)) then
                    bcgroupneeded = .true.
                    exit bcgroupfamloop
                end if
            end do bcgroupfamloop
 
            if (.not. bcgroupneeded) then
                cycle
            end if
 
            ! Now we know we *actually* need something from this BCGroup.
 
            ! Loop over each dimension individually since we have a scalar
            ! scatter.
            dimloop: do idim = 1, 3
 
                call vecgetarrayf90(exch%nodeValLocal, localptr, ierr)
                call echk(ierr, __file__, __line__)
 
                ! The local Pointer is just the localRandSurface we've set
                ! above.
                do j = 1, size(localptr)
                    localptr(i) = randsurface(3 * (j - 1) + idim)
                end do
 
                ! Restore the pointer
                call vecrestorearrayf90(exch%nodeValLocal, localptr, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Now scatter this to the zipper
                call vecscatterbegin(zipper%scatter, exch%nodeValLocal, &
                                     zipper%localVal, insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                call vecscatterend(zipper%scatter, exch%nodeValLocal, &
                                   zipper%localVal, insert_values, scatter_forward, ierr)
                call echk(ierr, __file__, __line__)
 
                ! The values we need are precisely what is in zipper%localVal
                call vecgetarrayf90(zipper%localVal, localptr, ierr)
                call echk(ierr, __file__, __line__)
 
                ! Just copy the received seeds into the random aray
                do j = 1, size(localptr)
                    ! Careful here becuase we have to interlate the dim
                    randsurface(3 * ii + 3 * (j - 1) + idim) = localptr(j)
                end do
 
            end do dimloop
 
            ! Increcment the running ii counter.
            ii = ii + size(localptr)
        end do
    end subroutine getsurfaceperturbation
 
end module warping