adjointUtils.F90

Go to the documentation of this file.

! This is a special function that is sued to alloc derivative values
! in blockpointers_d for use with the AD code.
module adjointutils
 
contains
 
    subroutine setupstateresidualmatrix(matrix, useAD, usePC, useTranspose, &
                                        useObjective, frozenTurb, level, useTurbOnly, useCoarseMats)
 
        !      Compute the state derivative matrix using a forward mode calc
        !      There are three different flags that determine how this
        !      routine is run:
        !      useAD: if True, AD is used for derivative calculation, if
        !             False, FD is used.
        !      usePC: if True, the reduced 1st order stencil with dissipation
        !             lumping is assembled instead of the actual exact
        !             full stencil jacobian
        !      useTranspose: If true, the transpose of dRdw is assembled.
        !                    For use with the adjoint this must be true.
        !      useObjective: If true, the force matrix is assembled
        !      level : What level to use to form the matrix. Level 1 is
        !              always the finest level
        !
        use block, only: flowdomsd
        use blockpointers
        use inputdiscretization
        use inputtimespectral
        use inputphysics
        use iteration
        use flowvarrefstate
        use inputadjoint
        use stencils
        use diffsizes
        use communication
        use adjointvars
        use turbmod
        use surfacefamilies, only: fullfamlist
        use oversetutilities, only: fractoweights
        use utils, only: echk, setpointers, getdirangle, setpointers_d
        use haloexchange, only: whalo2
        use masterroutines, only: block_res_state, master
        use amg, only: amglevels, a, coarseindices, coarseoversetindices
#ifndef USE_COMPLEX
        use masterroutines, only: block_res_state_d
#endif
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! PETSc Matrix Variable
        mat :: matrix
 
        ! Input Variables
        logical, intent(in) :: useAD, usePC, useTranspose, useObjective, frozenTurb
        logical, intent(in), optional :: useTurbOnly, useCoarseMats
        integer(kind=intType), intent(in) :: level
 
        ! Local variables.
        integer(kind=intType) :: ierr, nn, sps, sps2, i, j, k, l, ll, ii, jj, kk
        integer(kind=intType) :: nColor, iColor, jColor, irow, icol, fmDim, frow
        integer(kind=intType) :: nTransfer, nState, lStart, lEnd, tmp, icount, cols(8), nCol
        integer(kind=intType) :: n_stencil, i_stencil, m, iFringe, fInd, lvl, orderturbsave
        integer(kind=intType), dimension(:, :), pointer :: stencil
        real(kind=alwaysrealtype) :: delta_x, one_over_dx
        real(kind=realtype) :: weights(8), acousticscalesave
        real(kind=realtype), dimension(:, :), allocatable :: blk
        integer(kind=intType), dimension(2:10) :: coarseRows
        integer(kind=intType), dimension(8, 2:10) :: coarseCols
        integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, mm, colInd
        logical :: resetToRANS, turbOnly, flowRes, turbRes, buildCoarseMats
 
        ! Determine if we are assembling a turb only PC
        turbonly = .false.
        if (present(useturbonly)) then
            turbonly = useturbonly
        end if
 
        buildcoarsemats = .false.
        if (present(usecoarsemats)) then
            buildcoarsemats = usecoarsemats
        end if
 
        if (turbonly) then
            ! we are making a PC for turbulence only KSP
            flowres = .false.
            turbres = .true.
            lstart = nt1
            lend = nt2
            nstate = nt2 - nt1 + 1
        else
            ! We are making a "matrix" for either NK or adjoint
            ! Setup number of state variable based on turbulence assumption
            if (frozenturb) then
                flowres = .true.
                turbres = .false.
                lstart = 1
                lend = nwf
                nstate = nwf
            else
                flowres = .true.
                turbres = .true.
                lstart = 1
                lend = nw
                nstate = nw
            end if
        end if
 
        ! Generic block to use while setting values
        allocate (blk(nstate, nstate))
 
        ! Exchange data and call the residual to make sure its up to date
        ! withe current w
        call whalo2(1_inttype, 1_inttype, nw, .true., .true., .true.)
 
        ! This routine will not use the extra variables to block_res or the
        ! extra outputs, so we must zero them here
        alphad = zero
        betad = zero
        machd = zero
        machgridd = zero
        machcoefd = zero
        pointrefd = zero
        lengthrefd = zero
        pinfdimd = zero
        rhoinfdimd = zero
        tinfdimd = zero
 
        rkstage = 0
 
        ! Zero out the matrix before we start
        call matzeroentries(matrix, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the diagonal to 1 if the cell is not a compute cell:
 
        ! Make an identity block
        blk = zero
        do i = 1, nstate
            blk(i, i) = one
        end do
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            if (iblank(i, j, k) /= 1) then
                                irow = flowdoms(nn, level, sps)%globalCell(i, j, k)
                                cols(1) = irow
                                ncol = 1
 
                                if (buildcoarsemats) then
                                    do lvl = 1, amglevels - 1
                                        coarserows(lvl + 1) = coarseindices(nn, lvl)%arr(i, j, k)
                                        coarsecols(1, lvl + 1) = coarserows(lvl + 1)
                                    end do
                                end if
 
                                call setblock(blk)
                            end if
                        end do
                    end do
                end do
            end do
        end do
 
        ! Set a pointer to the correct set of stencil depending on if we are
        ! using the first order stencil or the full jacobian
 
        if (usepc) then
            if (viscous .and. viscpc) then
                stencil => visc_pc_stencil
                n_stencil = n_visc_pc
            else
                stencil => euler_pc_stencil
                n_stencil = n_euler_pc
            end if
 
            ! Very important to use only Second-Order dissipation for PC
            lumpeddiss = .true.
            ! We also do not apply acoustic scaling because artificial dissipation stabilizes the ILU factorization.
            ! This is mentioned in "Newton-Krylov-Schwarz Methods for Aerodynamics Problems: Compressible
            ! and Incompressible Flows on Unstructured Grids" by D. K. Kaushik, D. E. Keyes, and B. F. Smith (1998).
            ! The linear system will not converge if we reduce the artificial dissipation for the PC.
            acousticscalesave = acousticscalefactor
            acousticscalefactor = one
            ! also use first order advection terms for turbulence
            orderturbsave = orderturb
            orderturb = firstorder
        else
            if (viscous) then
                stencil => visc_drdw_stencil
                n_stencil = n_visc_drdw
            else
                stencil => euler_drdw_stencil
                n_stencil = n_euler_drdw
            end if
        end if
 
        ! Need to trick the residual evalution to use coupled (mean flow and
        ! turbulent) together.
 
        ! If we want to do the matrix on a coarser level, we must first
        ! restrict the fine grid solutions, since it is possible the
        ! NKsolver was used an the coarse grid solutions are (very!) out of
        ! date.
 
        ! Assembling matrix on coarser levels is not entirely implemented yet.
        currentlevel = level
        groundlevel = level
 
        ! Set delta_x
        delta_x = 1e-9_realtype
        one_over_dx = one / delta_x
        rkstage = 0
 
        ! Determine if we want to use frozenTurbulent Adjoint
        resettorans = .false.
        if (frozenturb .and. equations == ransequations) then
            equations = nsequations
            resettorans = .true.
        end if
 
        ! Allocate the additional memory we need for doing forward mode AD
        !  derivatives and copy any required reference values:
        if (.not. derivvarsallocated .and. usead) then
            call allocderivativevalues(level)
        end if
 
        ! For AD the initial seeds must all be zeroed.
        if (usead) then
            do nn = 1, ndom
                do sps = 1, ntimeintervalsspectral
                    call setpointers(nn, level, sps)
                    call zeroadseeds(nn, level, sps)
                end do
            end do
        end if
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
 
                ! Allocate some extra routines used only for assembly
                allocate ( &
                    flowdoms(nn, level, sps)%dw_deriv(2:il, 2:jl, 2:kl, 1:nw, 1:nw), &
                    flowdoms(nn, level, sps)%wtmp(0:ib, 0:jb, 0:kb, 1:nw), &
                    flowdoms(nn, level, sps)%dwtmp(0:ib, 0:jb, 0:kb, 1:nw), &
                    flowdoms(nn, level, sps)%dwtmp2(0:ib, 0:jb, 0:kb, 1:nw), &
                    stat=ierr)
                call echk(ierr, __file__, __line__)
 
                if (sps == 1) then
                    allocate (flowdoms(nn, level, sps)%color(0:ib, 0:jb, 0:kb), stat=ierr)
                    call echk(ierr, __file__, __line__)
                end if
            end do
        end do
 
        ! For the PC we don't linearize the shock sensor so it must be
        ! computed here.
 
        if (usepc) then
            call referenceshocksensor
        end if
 
        ! For FD, the initial reference values must be computed and stored.
        if (.not. usead) then
            call setfdreference(level)
        end if
 
        ! Master Domain Loop
        domainloopad: do nn = 1, ndom
 
            ! Set pointers to the first timeInstance...just to getSizes
            call setpointers(nn, level, 1)
            ! Set unknown sizes in diffSizes for AD routine
            isize1ofdrfbcdata = nbocos
            isize1ofdrfviscsubface = nviscbocos
 
            ! Setup the coloring for this block depending on if its
            ! drdw or a PC
 
            ! List of all Coloring Routines:
            !   Debugging Colorings Below:
            !       call setup_3x3x3_coloring(nn, level,  nColor)
            !       call setup_5x5x5_coloring(nn, level,  nColor)
            !       call setup_BF_coloring(nn, level,  nColor)
            !   Regular:
            !       call setup_PC_coloring(nn, level,  nColor)
            !       call setup_dRdw_euler_coloring(nn, level,  nColor)
            !       call setup_dRdw_visc_coloring(nn, level,  nColor)
 
            if (usepc) then
                if (viscous .and. viscpc) then
                    call setup_3x3x3_coloring(nn, level, ncolor) ! dense 3x3x3 coloring
                else
                    call setup_pc_coloring(nn, level, ncolor) ! Euler Colorings
                end if
            else
                if (viscous) then
                    !call setup_5x5x5_coloring(nn, level,  nColor)
                    call setup_drdw_visc_coloring(nn, level, ncolor)! Viscous/RANS
                else
                    call setup_drdw_euler_coloring(nn, level, ncolor) ! Euler Colorings
                end if
            end if
 
            spectralloop: do sps = 1, ntimeintervalsspectral
                ! Set pointers and (possibly derivative pointers)
                if (usead) then
                    call setpointers_d(nn, level, sps)
                else
                    call setpointers(nn, level, sps)
                end if
 
                ! Do Coloring and perturb states
                colorloop: do icolor = 1, ncolor
                    do sps2 = 1, ntimeintervalsspectral
                        flowdoms(nn, 1, sps2)%dw_deriv(:, :, :, :, :) = zero
                    end do
 
                    ! Master State Loop
                    stateloop: do l = lstart, lend
 
                        ! Reset All States and possibe AD seeds
                        do sps2 = 1, ntimeintervalsspectral
                            if (.not. usead) then
                                do ll = 1, nw
                                    do k = 0, kb
                                        do j = 0, jb
                                            do i = 0, ib
                                                flowdoms(nn, level, sps2)%w(i, j, k, ll) = &
                                                    flowdoms(nn, 1, sps2)%wtmp(i, j, k, ll)
                                            end do
                                        end do
                                    end do
                                end do
                            end if
 
                            if (usead) then
                                flowdomsd(nn, 1, sps2)%w = zero ! This is actually w seed
                            end if
                        end do
 
                        ! Peturb w or set AD Seed according to iColor. Note:
                        ! Do NOT try to putt he useAD if check inside the
                        ! color if check. ifort barfs hard-core on that and it
                        ! segfaults with AVX2
                        if (usead) then
                            do k = 0, kb
                                do j = 0, jb
                                    do i = 0, ib
                                        if (flowdoms(nn, 1, 1)%color(i, j, k) == icolor) then
                                            flowdomsd(nn, 1, sps)%w(i, j, k, l) = one
                                        end if
                                    end do
                                end do
                            end do
                        else
                            do k = 0, kb
                                do j = 0, jb
                                    do i = 0, ib
                                        if (flowdoms(nn, 1, 1)%color(i, j, k) == icolor) then
                                            w(i, j, k, l) = w(i, j, k, l) + delta_x
                                        end if
                                    end do
                                end do
                            end do
                        end if
 
                        ! Run Block-based residual
                        if (usead) then
#ifndef USE_COMPLEX
                            call block_res_state_d(nn, sps)
#else
                            print *, 'Forward AD routines are not complexified'
                            stop
#endif
                        else
                            call block_res_state(nn, sps, useflowres=flowres, useturbres=turbres)
                        end if
 
                        ! Set the computed residual in dw_deriv. If using FD,
                        ! actually do the FD calculation if AD, just copy out dw
                        ! in flowdomsd
 
                        ! Compute/Copy all derivatives
                        do sps2 = 1, ntimeintervalsspectral
                            do ll = lstart, lend
                                do k = 2, kl
                                    do j = 2, jl
                                        do i = 2, il
                                            if (usead) then
                                                flowdoms(nn, 1, sps2)%dw_deriv(i, j, k, ll, l) = &
                                                    flowdomsd(nn, 1, sps2)%dw(i, j, k, ll)
                                            else
                                                if (sps2 == sps) then
                                                    ! If the peturbation is on this
                                                    ! instance, we've computed the spatial
                                                    ! contribution so subtrace dwtmp
 
                                                    flowdoms(nn, 1, sps2)%dw_deriv(i, j, k, ll, l) = &
                                                        one_over_dx * &
                                                        (flowdoms(nn, 1, sps2)%dw(i, j, k, ll) - &
                                                         flowdoms(nn, 1, sps2)%dwtmp(i, j, k, ll))
                                                else
                                                    ! If the peturbation is on an off
                                                    ! instance, only subtract dwtmp2
                                                    ! which is the reference result
                                                    ! after initres
 
                                                    flowdoms(nn, 1, sps2)%dw_deriv(i, j, k, ll, l) = &
                                                        one_over_dx * ( &
                                                        flowdoms(nn, 1, sps2)%dw(i, j, k, ll) - &
                                                        flowdoms(nn, 1, sps2)%dwtmp2(i, j, k, ll))
                                                end if
                                            end if
                                        end do
                                    end do
                                end do
                            end do
                        end do
                    end do stateloop
 
                    ! Set derivatives by block in "matrix" after we've peturbed
                    ! all states in "color"
 
                    kloop: do k = 0, kb
                        jloop: do j = 0, jb
                            iloop: do i = 0, ib
                                colblank: if (flowdoms(nn, level, sps)%iblank(i, j, k) == 1 .or. &
                                              flowdoms(nn, level, sps)%iBlank(i, j, k) == -1) then
 
                                    ! If the cell we perturned ('iCol') is an
                                    ! interpolated cell, we don't actually use
                                    ! iCol, rather we use the 8 real donors that
                                    ! comprise the cell's value.
                                    if (flowdoms(nn, level, sps)%iblank(i, j, k) == 1) then
                                        cols(1) = flowdoms(nn, level, sps)%globalCell(i, j, k)
                                        ncol = 1
 
                                        if (buildcoarsemats) then
                                            do lvl = 1, amglevels - 1
                                                coarsecols(1, lvl + 1) = coarseindices(nn, lvl)%arr(i, j, k)
                                            end do
                                        end if
 
                                    else
                                        do m = 1, 8
                                            cols(m) = flowdoms(nn, level, sps)%gInd(m, i, j, k)
 
                                            if (buildcoarsemats) then
                                                do lvl = 1, amglevels - 1
                                                    coarsecols(m, lvl + 1) = &
                                                        coarseoversetindices(nn, lvl)%arr(m, i, j, k)
                                                end do
                                            end if
                                        end do
 
                                        find = fringeptr(1, i, j, k)
                                        call fractoweights(flowdoms(nn, level, sps)%fringes(find)%donorFrac, &
                                                           weights)
                                        ncol = 8
                                    end if
 
                                    colorcheck: if (flowdoms(nn, 1, 1)%color(i, j, k) == icolor) then
 
                                        ! i, j, k are now the "Center" cell that we
                                        ! actually petrubed. From knowledge of the
                                        ! stencil, we can simply take this cell and
                                        ! using the stencil, set the values around it
                                        ! in PETSc
 
                                        stencilloop: do i_stencil = 1, n_stencil
                                            ii = stencil(i_stencil, 1)
                                            jj = stencil(i_stencil, 2)
                                            kk = stencil(i_stencil, 3)
 
                                            ! Check to see if the cell in this
                                            ! sentcil is on a physical cell, not a
                                            ! halo/BC halo
                                            onblock: if (i + ii >= 2 .and. i + ii <= il .and. &
                                                         j + jj >= 2 .and. j + jj <= jl .and. &
                                                         k + kk >= 2 .and. k + kk <= kl) then
 
                                                irow = flowdoms(nn, level, sps)%globalCell( &
                                                       i + ii, j + jj, k + kk)
 
                                                if (buildcoarsemats) then
                                                    do lvl = 1, amglevels - 1
                                                        coarserows(lvl + 1) = &
                                                            coarseindices(nn, lvl)%arr(i + ii, j + jj, k + kk)
                                                    end do
                                                end if
 
                                                rowblank: if (flowdoms(nn, level, sps)% &
                                                              iblank(i + ii, j + jj, k + kk) == 1) then
 
                                                    centercell: if (ii == 0 .and. jj == 0 .and. kk == 0) then
                                                        usediagpc: if (usepc .and. usediagtspc) then
                                                            ! If we're doing the PC and we want
                                                            ! to use TS diagonal form, only set
                                                            ! values for on-time insintance
                                                            blk = flowdoms(nn, 1, sps)%dw_deriv(i + ii, &
                                                                                                j + jj, k + kk, &
                                                                                                lstart:lend, &
                                                                                                lstart:lend)
                                                            call setblock(blk)
                                                        else
                                                            ! Otherwise loop over spectral
                                                            ! instances and set all.
                                                            do sps2 = 1, ntimeintervalsspectral
                                                                irow = flowdoms(nn, level, sps2)% &
                                                                       globalcell(i + ii, &
                                                                                  j + jj, k + kk)
                                                                blk = flowdoms(nn, 1, sps2)% &
                                                                      dw_deriv(i + ii, &
                                                                               j + jj, k + kk, &
                                                                               lstart:lend, lstart:lend)
                                                                call setblock(blk)
                                                            end do
                                                        end if usediagpc
                                                    else
                                                        ! ALl other cells just set.
                                                        blk = flowdoms(nn, 1, sps)%dw_deriv(i + ii, j + jj, k + kk, &
                                                                                            lstart:lend, lstart:lend)
                                                        call setblock(blk)
                                                    end if centercell
                                                end if rowblank
                                            end if onblock
                                        end do stencilloop
                                    end if colorcheck
                                end if colblank
                            end do iloop
                        end do jloop
                    end do kloop
                end do colorloop
            end do spectralloop
        end do domainloopad
 
        ! PETSc Matrix Assembly begin
        call matassemblybegin(matrix, mat_final_assembly, ierr)
        call echk(ierr, __file__, __line__)
 
        if (buildcoarsemats) then
            do lvl = 2, amglevels
                call matassemblybegin(a(lvl), mat_final_assembly, ierr)
                call echk(ierr, __file__, __line__)
            end do
        end if
 
        ! Maybe we can do something useful while the communication happens?
        ! Deallocate the temporary memory used in this routine.
 
        ! Deallocate and reset values
        if (.not. usead) then
            call resetfdreference(level)
        end if
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                deallocate ( &
                    flowdoms(nn, 1, sps)%dw_deriv, &
                    flowdoms(nn, 1, sps)%wTmp, &
                    flowdoms(nn, 1, sps)%dwTmp, &
                    flowdoms(nn, 1, sps)%dwTmp2)
                if (sps == 1) then
                    deallocate (flowdoms(nn, 1, 1)%color)
                end if
            end do
        end do
 
        ! Return dissipation Parameters to normal -> VERY VERY IMPORTANT
        if (usepc) then
            lumpeddiss = .false.
            acousticscalefactor = acousticscalesave
            ! also recover the turbulence advection order
            orderturb = orderturbsave
        end if
 
        ! Reset the correct equation parameters if we were useing the frozen
        ! Turbulent
        if (resettorans) then
            equations = ransequations
        end if
 
        deallocate (blk)
 
        ! Complete the matrix assembly.
        call matassemblyend(matrix, mat_final_assembly, ierr)
        call echk(ierr, __file__, __line__)
 
        if (buildcoarsemats) then
            do lvl = 2, amglevels
                call matassemblyend(a(lvl), mat_final_assembly, ierr)
                call echk(ierr, __file__, __line__)
            end do
        end if
 
        call matsetoption(matrix, mat_new_nonzero_locations, petsc_false, ierr)
        call echk(ierr, __file__, __line__)
 
        ! ================= Important ===================
 
        ! We must run the residual computation to make sure that all
        ! intermediate variables are up to date. We can just call master
        ! for this. No need to recompute spatial terms.
        call master(.false.)
 
    contains
 
        subroutine setblock(blk)
            ! Sets a block at irow, icol, if useTranspose is False
            ! Sets a block at icol, irow with transpose of blk if useTranspose is True
 
            use genericisnan, only: myisnan
            implicit none
            real(kind=realtype), dimension(nState, nState) :: blk
 
            ! local variables
            integer(kind=intType) :: i, j, tmp, iRowSet, iColSet
            logical :: zeroFlag
            zeroflag = .false.
 
#ifndef USE_COMPLEX
            ! Check if the blk is all zeros
 
            zeroflag = .true.
            do i = 1, nstate
                do j = 1, nstate
                    if (.not. blk(i, j) == zero) then
                        zeroflag = .false.
                    end if
                end do
            end do
 
            ! Check if the blk has nan
            if (myisnan(sum(blk))) then
                print *, 'Bad Block:', blk
                print *, 'irow:', irow
                print *, 'icol', cols(1:ncol)
                print *, 'nn:', nn
                print *, 'ijk:', i, j, k
                call echk(1, __file__, __line__)
            end if
#endif
 
            if (.not. zeroflag) then
                if (ncol == 1) then
                    if (usetranspose) then
                        blk = transpose(blk)
                        call matsetvaluesblocked(matrix, 1, cols(1), 1, irow, blk, &
                                                 add_values, ierr)
                        call echk(ierr, __file__, __line__)
                    else
                        call matsetvaluesblocked(matrix, 1, irow, 1, cols(1), blk, &
                                                 add_values, ierr)
                        call echk(ierr, __file__, __line__)
                    end if
                else
                    if (usetranspose) then
                        blk = transpose(blk)
                        do m = 1, ncol
                            if (cols(m) >= 0) then
                                call matsetvaluesblocked(matrix, 1, cols(m), 1, irow, blk * weights(m), &
                                                         add_values, ierr)
                                call echk(ierr, __file__, __line__)
                            end if
                        end do
                    else
                        do m = 1, ncol
                            if (cols(m) >= 0) then
                                call matsetvaluesblocked(matrix, 1, irow, 1, cols(m), blk * weights(m), &
                                                         add_values, ierr)
                                call echk(ierr, __file__, __line__)
                            end if
                        end do
                    end if
                end if
 
                ! Extension for setting coarse grids:
                if (buildcoarsemats) then
                    if (ncol == 1) then
                        do lvl = 2, amglevels
                            if (usetranspose) then
                                ! Loop over the coarser levels
                                call matsetvaluesblocked(a(lvl), 1, coarsecols(1, lvl), 1, coarserows(lvl), &
                                                         blk, add_values, ierr)
                            else
                                call matsetvaluesblocked(a(lvl), 1, coarserows(lvl), 1, coarsecols(1, lvl), &
                                                         blk, add_values, ierr)
                            end if
                        end do
                    else
                        do m = 1, ncol
                            do lvl = 2, amglevels
                                if (coarsecols(m, lvl) >= 0) then
                                    if (usetranspose) then
                                        ! Loop over the coarser levels
                                        call matsetvaluesblocked(a(lvl), 1, coarsecols(m, lvl), 1, coarserows(lvl), &
                                                                 blk * weights(m), add_values, ierr)
                                    else
                                        call matsetvaluesblocked(a(lvl), 1, coarserows(lvl), 1, coarsecols(m, lvl), &
                                                                 blk * weights(m), add_values, ierr)
                                    end if
                                end if
                            end do
                        end do
                    end if
                end if
            end if
        end subroutine setblock
    end subroutine setupstateresidualmatrix
 
    subroutine allocderivativevalues(level)
 
        use constants
        use block, only: flowdoms, flowdomsd
        use blockpointers, only: ndom, il, jl, kl, ie, je, ke, ib, jb, kb, bcdata, &
                                 nbocos, nviscbocos
        use inputtimespectral, only: ntimeintervalsspectral
        use flowvarrefstate, only: winf, winfd, nw, nt1, nt2
        use inputdiscretization, only: useapproxwalldistance
        use inputphysics, only: walldistanceneeded
        use communication, only: adflow_comm_world
        use walldistancedata, only: xsurfvec, xsurfvecd!, PETSC_DETERMINE
        use bcpointers_b
        use adjointvars, only: derivvarsallocated
        use utils, only: echk, setpointers, getdirangle
        use cgnsgrid, only: cgnsdoms, cgnsdomsd, cgnsndom
        implicit none
 
        ! Input parameters
        integer(kind=intType) :: level
 
        ! Local variables
        integer(kind=intType) :: sps, ierr, i, j, k, l, mm, nn
        integer(kind=intType) :: iBeg, jBeg, iStop, jStop, isizemax, jsizemax
        integer(kind=intType) :: inBeg, jnBeg, inStop, jnStop
        integer(kind=intType) :: massShape(2), max_face_size
        integer(kind=intType) :: iBoco, nDataset, iData, nDirichlet, iDirichlet, nArray
        ! First create the derivative flowdoms structure flowDomsd. Note we
        ! only allocate information for the finest grid.
 
        allocate (flowdomsd(ndom, 1, ntimeintervalsspectral), stat=ierr)
        call echk(ierr, __file__, __line__)
 
        ! winfd hasn't be allocated so we'll do it here
        allocate (winfd(size(winf)))
 
        ! If we are not using RANS with walDistance create a dummy xSurfVec
        ! since one does not yet exist
        if (.not. walldistanceneeded .or. .not. useapproxwalldistance) then
            do sps = 1, ntimeintervalsspectral
                call veccreatempi(adflow_comm_world, 1, petsc_determine, xsurfvec(1, sps), ierr)
                call echk(ierr, __file__, __line__)
            end do
        end if
 
        ! Duplicate the PETSc Xsurf Vec, but only on the first level:
        allocate (xsurfvecd(ntimeintervalsspectral))
        do sps = 1, ntimeintervalsspectral
            call vecduplicate(xsurfvec(1, sps), xsurfvecd(sps), ierr)
            call echk(ierr, __file__, __line__)
        end do
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
 
                ! Allocate d2wall if not already done so
                if (.not. associated(flowdoms(nn, 1, sps)%d2wall)) then
                    allocate (flowdoms(nn, 1, sps)%d2wall(2:il, 2:jl, 2:kl))
                    call echk(ierr, __file__, __line__)
                end if
 
                ! Now allocate all valus that have a differentiable
                ! dependence.
                allocate ( &
                    flowdomsd(nn, level, sps)%x(0:ie, 0:je, 0:ke, 3), &
                    flowdomsd(nn, level, sps)%vol(0:ib, 0:jb, 0:kb), &
                    flowdomsd(nn, level, sps)%si(0:ie, 1:je, 1:ke, 3), &
                    flowdomsd(nn, level, sps)%sj(1:ie, 0:je, 1:ke, 3), &
                    flowdomsd(nn, level, sps)%sk(1:ie, 1:je, 0:ke, 3), &
                    flowdomsd(nn, level, sps)%rotMatrixI(il, 2:jl, 2:kl, 3, 3), &
                    flowdomsd(nn, level, sps)%rotMatrixJ(2:il, jl, 2:kl, 3, 3), &
                    flowdomsd(nn, level, sps)%rotMatrixK(2:il, 2:jl, kl, 3, 3), &
                    flowdomsd(nn, level, sps)%s(ie, je, ke, 3), &
                    flowdomsd(nn, level, sps)%sFaceI(0:ie, je, ke), &
                    flowdomsd(nn, level, sps)%sFaceJ(ie, 0:je, ke), &
                    flowdomsd(nn, level, sps)%sFaceK(ie, je, 0:ke), &
                    flowdomsd(nn, level, sps)%w(0:ib, 0:jb, 0:kb, 1:nw), &
                    flowdomsd(nn, level, sps)%dw(0:ib, 0:jb, 0:kb, 1:nw), &
                    flowdomsd(nn, level, sps)%fw(0:ib, 0:jb, 0:kb, 1:nw), &
                    flowdomsd(nn, level, sps)%scratch(0:ib, 0:jb, 0:kb, 5), &
                    flowdomsd(nn, level, sps)%p(0:ib, 0:jb, 0:kb), &
                    flowdomsd(nn, level, sps)%gamma(0:ib, 0:jb, 0:kb), &
                    flowdomsd(nn, level, sps)%aa(0:ib, 0:jb, 0:kb), &
                    flowdomsd(nn, level, sps)%ux(il, jl, kl), &
                    flowdomsd(nn, level, sps)%uy(il, jl, kl), &
                    flowdomsd(nn, level, sps)%uz(il, jl, kl), &
                    flowdomsd(nn, level, sps)%vx(il, jl, kl), &
                    flowdomsd(nn, level, sps)%vy(il, jl, kl), &
                    flowdomsd(nn, level, sps)%vz(il, jl, kl), &
                    flowdomsd(nn, level, sps)%wx(il, jl, kl), &
                    flowdomsd(nn, level, sps)%wy(il, jl, kl), &
                    flowdomsd(nn, level, sps)%wz(il, jl, kl), &
                    flowdomsd(nn, level, sps)%qx(il, jl, kl), &
                    flowdomsd(nn, level, sps)%qy(il, jl, kl), &
                    flowdomsd(nn, level, sps)%qz(il, jl, kl), &
                    flowdomsd(nn, level, sps)%rlv(0:ib, 0:jb, 0:kb), &
                    flowdomsd(nn, level, sps)%rev(0:ib, 0:jb, 0:kb), &
                    flowdomsd(nn, level, sps)%dtl(1:ie, 1:je, 1:ke), &
                    flowdomsd(nn, level, sps)%radI(1:ie, 1:je, 1:ke), &
                    flowdomsd(nn, level, sps)%radJ(1:ie, 1:je, 1:ke), &
                    flowdomsd(nn, level, sps)%radK(1:ie, 1:je, 1:ke), &
                    flowdomsd(nn, level, sps)%BCData(nbocos), &
                    flowdomsd(nn, level, sps)%bmti1(je, ke, nt1:nt2, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bmti2(je, ke, nt1:nt2, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bmtj1(ie, ke, nt1:nt2, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bmtj2(ie, ke, nt1:nt2, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bmtk1(ie, je, nt1:nt2, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bmtk2(ie, je, nt1:nt2, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bvti1(je, ke, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bvti2(je, ke, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bvtj1(ie, ke, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bvtj2(ie, ke, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bvtk1(ie, je, nt1:nt2), &
                    flowdomsd(nn, level, sps)%bvtk2(ie, je, nt1:nt2), &
                    flowdomsd(nn, level, sps)%d2Wall(2:il, 2:jl, 2:kl), &
                    stat=ierr)
                call echk(ierr, __file__, __line__)
 
                allocate (flowdomsd(nn, level, sps)%viscSubface(nviscbocos), &
                          stat=ierr)
                call echk(ierr, __file__, __line__)
 
                ! Set the number of bocos/viscbocs
                flowdomsd(nn, level, sps)%nBocos = flowdoms(nn, level, sps)%nbocos
                flowdomsd(nn, level, sps)%nViscBocos = flowdoms(nn, level, sps)%nViscBocos
                bocoloop: do mm = 1, nbocos
 
                    ! Store the cell range of the boundary subface
                    ! a bit easier.
 
                    ibeg = bcdata(mm)%icbeg; istop = bcdata(mm)%icend
                    jbeg = bcdata(mm)%jcbeg; jstop = bcdata(mm)%jcend
 
                    inbeg = bcdata(mm)%inBeg; instop = bcdata(mm)%inEnd
                    jnbeg = bcdata(mm)%jnBeg; jnstop = bcdata(mm)%jnEnd
                    allocate ( &
                        flowdomsd(nn, level, sps)%BCData(mm)%norm(ibeg:istop, jbeg:jstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%rface(ibeg:istop, jbeg:jstop), &
                        flowdomsd(nn, level, sps)%BCData(mm)%Fp(inbeg + 1:instop, jnbeg + 1:jnstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%Fv(inbeg + 1:instop, jnbeg + 1:jnstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%Tp(inbeg:instop, jnbeg:jnstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%Tv(inbeg:instop, jnbeg:jnstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%F(inbeg:instop, jnbeg:jnstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%T(inbeg:instop, jnbeg:jnstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%area(inbeg + 1:instop, jnbeg + 1:jnstop), &
                        flowdomsd(nn, level, sps)%BCData(mm)%uSlip(ibeg:istop, jbeg:jstop, 3), &
                        flowdomsd(nn, level, sps)%BCData(mm)%TNS_Wall(ibeg:istop, jbeg:jstop), &
                        flowdomsd(nn, level, sps)%BCData(mm)%ptInlet(ibeg:istop, jbeg:jstop), &
                        flowdomsd(nn, level, sps)%BCData(mm)%htInlet(ibeg:istop, jbeg:jstop), &
                        flowdomsd(nn, level, sps)%BCData(mm)%ttInlet(ibeg:istop, jbeg:jstop), &
                        flowdomsd(nn, level, sps)%BCData(mm)%turbInlet(ibeg:istop, jbeg:jstop, nt1:nt2), &
                        flowdomsd(nn, level, sps)%BCData(mm)%ps(ibeg:istop, jbeg:jstop), stat=ierr)
 
                    call echk(ierr, __file__, __line__)
                end do bocoloop
 
                viscbocoloop: do mm = 1, nviscbocos
 
                    ibeg = bcdata(mm)%inBeg + 1; istop = bcdata(mm)%inEnd
                    jbeg = bcdata(mm)%jnBeg + 1; jstop = bcdata(mm)%jnEnd
 
                    allocate ( &
                        flowdomsd(nn, level, sps)%viscSubface(mm)%tau(ibeg:istop, jbeg:jstop, 6), &
                        flowdomsd(nn, level, sps)%viscSubface(mm)%q(ibeg:istop, jbeg:jstop, 6), &
                        stat=ierr)
                    call echk(ierr, __file__, __line__)
                end do viscbocoloop
            end do
        end do
 
        ! Allocate the derivatives values for the CGNS data structure used
        ! to store boundary condition values
        do nn = 1, cgnsndom
            nbocos = cgnsdoms(nn)%nBocos
            allocate (cgnsdomsd(nn)%bocoInfo(nbocos))
            do iboco = 1, nbocos
                if (associated(cgnsdoms(nn)%bocoInfo(iboco)%dataSet)) then
                    ndataset = size(cgnsdoms(nn)%bocoInfo(iboco)%dataSet)
                    allocate (cgnsdomsd(nn)%bocoInfo(iboco)%dataSet(ndataset))
 
                    do idata = 1, ndataset
                        if (associated(cgnsdoms(nn)%bocoInfo(iboco)%dataSet(idata)%dirichletArrays)) then
                            ndirichlet = size(cgnsdoms(nn)%bocoInfo(iboco)%dataSet(idata)%dirichletArrays)
                            allocate (cgnsdomsd(nn)%bocoInfo(iboco)%dataSet(idata)%dirichletArrays(ndirichlet))
 
                            do idirichlet = 1, ndirichlet
                                narray = size(cgnsdoms(nn)%bocoInfo(iboco)%dataSet(idata) &
                                              %dirichletArrays(idirichlet)%dataArr)
                                allocate (cgnsdomsd(nn)%bocoInfo(iboco)% &
                                          dataset(idata)%dirichletArrays(idirichlet)%dataArr(narray))
                                cgnsdomsd(nn)%bocoInfo(iboco)%dataSet(idata)% &
                                    dirichletarrays(idirichlet)%dataArr(narray) = zero
                            end do
                        end if
                    end do
                end if
            end do
        end do
 
        derivvarsallocated = .true.
    end subroutine allocderivativevalues
 
    subroutine zeroadseeds(nn, level, sps)
 
        use constants
        use block, only: flowdomsd, flowdoms
        use blockpointers
        use inputtimespectral
        use flowvarrefstate
        use inputphysics
        use bcpointers_b
        use communication
        use oversetdata, only: oversetpresent
        use cgnsgrid, only: cgnsdoms, cgnsdomsd, cgnsndom
        use actuatorregiondata, only: nactuatorregions, actuatorregionsd
        implicit none
 
        ! Input parameters
        integer(kind=intType) :: nn, level, sps
 
        ! Working parameters
        integer(kind=intType) :: mm, i, iDom
        integer(kind=intType) :: iBoco, iData, iDirichlet
        flowdomsd(nn, level, sps)%d2wall = zero
        flowdomsd(nn, level, sps)%x = zero
        flowdomsd(nn, level, sps)%si = zero
        flowdomsd(nn, level, sps)%sj = zero
        flowdomsd(nn, level, sps)%sk = zero
        flowdomsd(nn, level, sps)%vol = zero
 
        flowdomsd(nn, level, sps)%s = zero
        flowdomsd(nn, level, sps)%sFaceI = zero
        flowdomsd(nn, level, sps)%sFaceJ = zero
        flowdomsd(nn, level, sps)%sFaceK = zero
 
        flowdomsd(nn, level, sps)%w = zero
        flowdomsd(nn, level, sps)%dw = zero
        flowdomsd(nn, level, sps)%fw = zero
        flowdomsd(nn, level, sps)%scratch = zero
 
        flowdomsd(nn, level, sps)%p = zero
        flowdomsd(nn, level, sps)%gamma = zero
        flowdomsd(nn, level, sps)%aa = zero
 
        flowdomsd(nn, level, sps)%rlv = zero
        flowdomsd(nn, level, sps)%rev = zero
 
        flowdomsd(nn, level, sps)%dtl = zero
 
        flowdomsd(nn, level, sps)%radI = zero
        flowdomsd(nn, level, sps)%radJ = zero
        flowdomsd(nn, level, sps)%radK = zero
 
        flowdomsd(nn, level, sps)%ux = zero
        flowdomsd(nn, level, sps)%uy = zero
        flowdomsd(nn, level, sps)%uz = zero
        flowdomsd(nn, level, sps)%vx = zero
        flowdomsd(nn, level, sps)%vy = zero
        flowdomsd(nn, level, sps)%vz = zero
        flowdomsd(nn, level, sps)%wx = zero
        flowdomsd(nn, level, sps)%wy = zero
        flowdomsd(nn, level, sps)%wz = zero
        flowdomsd(nn, level, sps)%qx = zero
        flowdomsd(nn, level, sps)%qy = zero
        flowdomsd(nn, level, sps)%qz = zero
 
        flowdomsd(nn, level, sps)%bmti1 = zero
        flowdomsd(nn, level, sps)%bmti2 = zero
        flowdomsd(nn, level, sps)%bmtj1 = zero
        flowdomsd(nn, level, sps)%bmtj2 = zero
        flowdomsd(nn, level, sps)%bmtk1 = zero
        flowdomsd(nn, level, sps)%bmtk2 = zero
        flowdomsd(nn, level, sps)%bvti1 = zero
        flowdomsd(nn, level, sps)%bvti2 = zero
        flowdomsd(nn, level, sps)%bvtj1 = zero
        flowdomsd(nn, level, sps)%bvtj2 = zero
        flowdomsd(nn, level, sps)%bvtk1 = zero
        flowdomsd(nn, level, sps)%bvtk2 = zero
 
        bocoloop: do mm = 1, flowdoms(nn, level, sps)%nBocos
            flowdomsd(nn, level, sps)%BCData(mm)%norm = zero
            flowdomsd(nn, level, sps)%bcData(mm)%rface = zero
            flowdomsd(nn, level, sps)%bcData(mm)%Fv = zero
            flowdomsd(nn, level, sps)%bcData(mm)%Fp = zero
            flowdomsd(nn, level, sps)%bcData(mm)%Tv = zero
            flowdomsd(nn, level, sps)%bcData(mm)%Tp = zero
            flowdomsd(nn, level, sps)%bcData(mm)%area = zero
            flowdomsd(nn, level, sps)%BCData(mm)%uSlip = zero
            flowdomsd(nn, level, sps)%BCData(mm)%TNS_Wall = zero
            flowdomsd(nn, level, sps)%BCData(mm)%ptInlet = zero
            flowdomsd(nn, level, sps)%BCData(mm)%htInlet = zero
            flowdomsd(nn, level, sps)%BCData(mm)%ttInlet = zero
            flowdomsd(nn, level, sps)%BCData(mm)%turbInlet = zero
            flowdomsd(nn, level, sps)%BCData(mm)%ps = zero
        end do bocoloop
 
        viscbocoloop: do mm = 1, flowdoms(nn, level, sps)%nViscBocos
            flowdomsd(nn, level, sps)%viscSubface(mm)%tau = zero
            flowdomsd(nn, level, sps)%viscSubface(mm)%q = zero
        end do viscbocoloop
 
        ! For overset, the weights may be active in the comm structure. We
        ! need to zero them before we can accumulate.
        if (oversetpresent) then
            ! Pointers to the overset comms to make it easier to read
            sends: do i = 1, commpatternoverset(level, sps)%nProcSend
                commpatternoverset(level, sps)%sendList(i)%interpd = zero
            end do sends
            internaloverset(level, sps)%donorInterpd = zero
        end if
 
        alphad = zero
        betad = zero
        machd = zero
        machgridd = zero
        machcoefd = zero
        pinfdimd = zero
        tinfdimd = zero
        rhoinfdimd = zero
        rgasdimd = zero
        pointrefd = zero
        prefd = zero
        rhorefd = zero
        trefd = zero
        murefd = zero
        urefd = zero
        hrefd = zero
        timerefd = zero
        pinfd = zero
        pinfcorrd = zero
        rhoinfd = zero
        uinfd = zero
        rgasd = zero
        muinfd = zero
        gammainfd = zero
        winfd = zero
        veldirfreestreamd = zero
        liftdirectiond = zero
        dragdirectiond = zero
 
        ! Zero all the reverse seeds in the dirichlet input arrays
        do idom = 1, cgnsndom
            do iboco = 1, cgnsdoms(idom)%nBocos
                if (associated(cgnsdoms(idom)%bocoInfo(iboco)%dataSet)) then
                    do idata = 1, size(cgnsdoms(idom)%bocoInfo(iboco)%dataSet)
                        if (associated(cgnsdoms(idom)%bocoInfo(iboco)%dataSet(idata)%dirichletArrays)) then
                            do idirichlet = 1, &
                                size(cgnsdoms(idom)%bocoInfo(iboco)%dataSet(idata)%dirichletArrays)
                                cgnsdomsd(idom)%bocoInfo(iboco)%dataSet(idata)%dirichletArrays(idirichlet)%dataArr(:) &
                                    = zero
                            end do
                        end if
                    end do
                end if
            end do
        end do
 
        ! And the reverse seeds in the actuator zones
        do i = 1, nactuatorregions
            actuatorregionsd(i)%force = zero
            actuatorregionsd(i)%torque = zero
            actuatorregionsd(i)%heat = zero
            actuatorregionsd(i)%volume = zero
        end do
 
    end subroutine zeroadseeds
    ! This is a special function that is sued to dealloc derivative values
    ! in blockpointers_d for use with the AD code.
 
    ! This routine setups "coloring" stencils for various sizes
 
    subroutine setup_pc_coloring(nn, level, nColor)
 
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb
        use utils, only: setpointers
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, level
 
        ! Output parameters
        integer(kind=intTYpe), intent(out) :: nColor
 
        ! Working
        integer(kind=intType) :: i, j, k
 
        call setpointers(nn, level, 1)
        !DIR$ NOVECTOR
        do k = 0, kb
            do j = 0, jb
                do i = 0, ib
                    ! Add the extra one for 1-based numbering (as opposed to zero-based)
                    flowdoms(nn, level, 1)%color(i, j, k) = &
                        mod(i + 5 * j + 4 * k, 7) + 1
                end do
            end do
        end do
 
        ncolor = 7
 
    end subroutine setup_pc_coloring
 
    subroutine setup_drdw_euler_coloring(nn, level, nColor)
 
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb
        use utils, only: setpointers
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, level
 
        ! Output parameters
        integer(kind=intTYpe), intent(out) :: nColor
 
        ! Working
        integer(kind=intType) :: i, j, k
 
        call setpointers(nn, level, 1)
        !DIR$ NOVECTOR
        do k = 0, kb
            do j = 0, jb
                do i = 0, ib
                    ! Add the extra one for 1-based numbering (as opposed to zero-based)
                    flowdoms(nn, level, 1)%color(i, j, k) = &
                        mod(i + 3 * j + 4 * k, 13) + 1
                end do
            end do
        end do
 
        ncolor = 13
 
    end subroutine setup_drdw_euler_coloring
 
    subroutine setup_drdw_visc_coloring(nn, level, nColor)
 
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb
        use utils, only: setpointers
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, level
 
        ! Output parameters
        integer(kind=intTYpe), intent(out) :: nColor
 
        ! Working
        integer(kind=intType) :: i, j, k
 
        call setpointers(nn, level, 1) ! Just to get the correct sizes
        !DIR$ NOVECTOR
        do k = 0, kb
            do j = 0, jb
                do i = 0, ib
                    ! Add the extra one for 1-based numbering (as opposed to zero-based)
                    flowdoms(nn, level, 1)%color(i, j, k) = &
                        mod(i + 19 * j + 11 * k, 35) + 1
                end do
            end do
        end do
 
        ncolor = 35
 
    end subroutine setup_drdw_visc_coloring
 
    ! -------------------------------------------------------------
    !                   Debugging Color Colorings
    ! -------------------------------------------------------------
 
    subroutine setup_3x3x3_coloring(nn, level, nColor)
 
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb
        use utils, only: setpointers
        implicit none
 
        ! This is a dense 3x3x3 cube for debugging only
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, level
 
        ! Output parameters
        integer(kind=intTYpe), intent(out) :: nColor
 
        ! Working
        integer(kind=intType) :: i, j, k, modi, modj, modk
 
        call setpointers(nn, level, 1)
        !DIR$ NOVECTOR
        do k = 0, kb
            do j = 0, jb
                do i = 0, ib
                    ! Add the extra one for 1-based numbering (as opposed to zero-based)
                    modi = mod(i, 3)
                    modj = mod(j, 3)
                    modk = mod(k, 3)
 
                    flowdoms(nn, level, 1)%color(i, j, k) = modi + 3 * modj + 9 * modk + 1
 
                end do
            end do
        end do
 
        ncolor = 27
    end subroutine setup_3x3x3_coloring
 
    subroutine setup_5x5x5_coloring(nn, level, nColor)
 
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb
        use utils, only: setpointers
        implicit none
 
        ! This is a dense 5x5x5 cube for debugging only
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, level
 
        ! Output parameters
        integer(kind=intTYpe), intent(out) :: nColor
 
        ! Working
        integer(kind=intType) :: i, j, k, modi, modj, modk
 
        call setpointers(nn, level, 1)
        !DIR$ NOVECTOR
        do k = 0, kb
            do j = 0, jb
                do i = 0, ib
                    ! Add the extra one for 1-based numbering (as opposed to zero-based)
                    modi = mod(i, 5)
                    modj = mod(j, 5)
                    modk = mod(k, 5)
 
                    flowdoms(nn, level, 1)%color(i, j, k) = modi + 5 * modj + 25 * modk + 1
 
                end do
            end do
        end do
 
        ncolor = 125
    end subroutine setup_5x5x5_coloring
 
    subroutine setup_bf_coloring(nn, level, nColor)
 
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb
        use utils, only: setpointers
        implicit none
 
        ! Input parameters
        integer(kind=intType), intent(in) :: nn, level
 
        ! Output parameters
        integer(kind=intTYpe), intent(out) :: nColor
 
        ! Working
        integer(kind=intType) :: i, j, k
 
        ! This is a REALLY brute force coloring for debugging
 
        call setpointers(nn, level, 1)
        !DIR$ NOVECTOR
        do k = 0, kb
            do j = 0, jb
                do i = 0, ib
                    ! Add the extra one for 1-based numbering (as opposed to zero-based)
 
                    flowdoms(nn, level, 1)%color(i, j, k) = i + j * (ib + 1) + k * ((ib + 1) * (jb + 1)) + 1
                end do
            end do
        end do
 
        ncolor = (ib + 1) * (jb + 1) * (kb + 1)
    end subroutine setup_bf_coloring
 
    subroutine mymatcreate(matrix, blockSize, m, n, nnzDiagonal, nnzOffDiag, &
                           file, line)
        ! Function to create petsc matrix to make stuff a little cleaner in
        ! the code above. Also, PETSc always thinks is a good idea to
        ! RANDOMLY change syntax between versions so this way there is only
        ! one place to make a change based on petsc version.
 
        use constants
        use communication, only: adflow_comm_world
        use utils, only: echk, setpointers
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        mat matrix
        integer(kind=intType), intent(in) :: blockSize, m, n
        integer(kind=intType), intent(in), dimension(*) :: nnzDiagonal, nnzOffDiag
        character(len=*) :: file
        integer(kind=intType) :: ierr, line
        ! if (blockSize > 1) then
        call matcreatebaij(adflow_comm_world, blocksize, &
                           m, n, petsc_determine, petsc_determine, &
                           0, nnzdiagonal, 0, nnzoffdiag, matrix, ierr)
        ! else
        ! call MatCreateAIJ(ADFLOW_COMM_WORLD,&
        ! m, n, PETSC_DETERMINE, PETSC_DETERMINE, &
        ! 0, nnzDiagonal, 0, nnzOffDiag, matrix, ierr)
        call echk(ierr, file, line)
        ! end if
 
        ! Warning: The array values is logically two-dimensional,
        ! containing the values that are to be inserted. By default the
        ! values are given in row major order, which is the opposite of
        ! the Fortran convention, meaning that the value to be put in row
        ! idxm[i] and column idxn[j] is located in values[i*n+j]. To allow
        ! the insertion of values in column major order, one can call the
        ! command MatSetOption(Mat A, MAT COLUMN ORIENTED);
 
        call matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
        call echk(ierr, __file__, __line__)
 
        call matsetoption(matrix, mat_new_nonzero_allocation_err, petsc_false, ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine mymatcreate
 
    subroutine mykspmonitor(myKsp, n, rnorm, dummy, ierr)
        !
        !      This is a user-defined routine for monitoring the KSP
        !      iterative solvers. Instead of outputing the L2-norm at every
        !      iteration (default PETSc monitor), it only does it every
        !      'adjMonStep' iterations.
        !
        use adjointpetsc
        use inputadjoint
        use communication
        implicit none
        !
        !     Subroutine arguments.
        !
        ! myKsp - Iterative context
        ! n     - Iteration number
        ! rnorm - 2-norm (preconditioned) residual value
        ! dummy - Optional user-defined monitor context (unused here)
        ! ierr  - Return error code
 
        real(kind=realtype), pointer, dimension(:, :) :: myksp
        integer(kind=intType) :: n, dummy, ierr
        real(kind=realtype) :: rnorm
 
        ! Write the residual norm to stdout every adjMonStep iterations.
 
        if (mod(n, adjmonstep) == 0) then
            if (myid == 0) write (*, "(I4, 1X, A, 1X, ES16.10)") n, 'KSP Residual norm', rnorm
        end if
 
        ierr = 0
 
    end subroutine mykspmonitor
 
    subroutine setupstandardksp(kspObject, kspObjectType, gmresRestart, preConSide, &
                                globalPCType, ASMOverlap, globalPreConIts, localPCType, &
                                localMatrixOrdering, localFillLevel, localPreConIts)
 
        ! This function sets up the supplied kspObject in the following
        ! specific fashion. The reason this setup is in
        ! its own function is that it is used in the following places:
        ! 1. Setting up the preconditioner to use for the ANK and NK solvers
        ! 2. Setting up the preconditioner to use for the adjoint solver
        ! 3. Setting up the smoothers on the coarse levels for algebraic multigrid
        !
        ! The hierarchy of the setup is:
        !  kspObject --> Supplied KSP object
        !  |
        !  --> master_PC --> Preconditioner type set to KSP
        !      |
        !      --> master_PC_KSP --> KSP type set to Richardson with 'globalPreConIts'
        !          |
        !           --> globalPC --> PC type set to 'globalPCType'
        !               |            Usually Additive Schwarz and overlap is set
        !               |            with 'ASMOverlap'. Use 0 to get BlockJacobi
        !               |
        !               --> subKSP --> KSP type set to Richardon with 'LocalPreConIts'
        !                   |
        !                   --> subPC -->  PC type set to 'loclaPCType'.
        !                                  Usually ILU. 'localFillLevel' is
        !                                  set and 'localMatrixOrder' is used.
        !
        ! Note that if globalPreConIts=1 then maser_PC_KSP is NOT created and master_PC=globalPC
        ! and if localPreConIts=1 then subKSP is set to preOnly.
        use constants
        use utils, only: echk
        use inputadjoint, only: gmresorthogtype
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Input Params
        ksp kspobject
        character(len=*), intent(in) :: kspObjectType, preConSide
        character(len=*), intent(in) :: globalPCType, localPCType
        character(len=*), intent(in) :: localMatrixOrdering
        integer(kind=intType), intent(in) :: ASMOverlap, localFillLevel, gmresRestart
        integer(kind=intType), intent(in) :: globalPreConIts, localPreConIts
 
        ! Working Variables
        pc master_pc, globalpc, subpc
        ksp master_pc_ksp, subksp
        integer(kind=intType) :: nlocal, first, ierr
 
        ! First, KSPSetFromOptions MUST be called
        call kspsetfromoptions(kspobject, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the type of solver to use:
        call kspsettype(kspobject, kspobjecttype, ierr)
        call echk(ierr, __file__, __line__)
 
        ! If we're using GMRES set the possible gmres restart
        call kspgmressetrestart(kspobject, gmresrestart, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the orthogonalization method for GMRES
        select case (gmresorthogtype)
        case ('modified_gram_schmidt')
            ! Use modified Gram-Schmidt
            call kspgmressetorthogonalization(kspobject, kspgmresmodifiedgramschmidtorthogonalization, ierr)
        case ('cgs_never_refine')
            ! Use classical Gram-Schmidt with no refinement
            call kspgmressetcgsrefinementtype(kspobject, ksp_gmres_cgs_refine_never, ierr)
        case ('cgs_refine_if_needed')
            ! Use classical Gram-Schmidt with refinement if needed
            call kspgmressetcgsrefinementtype(kspobject, ksp_gmres_cgs_refine_ifneeded, ierr)
        case ('cgs_always_refine')
            ! Use classical Gram-Schmidt with refinement at every iteration
            call kspgmressetcgsrefinementtype(kspobject, ksp_gmres_cgs_refine_always, ierr)
        end select
        call echk(ierr, __file__, __line__)
 
        ! Set the preconditioner side from option:
        if (trim(preconside) == 'right') then
            call kspsetpcside(kspobject, pc_right, ierr)
        else
            call kspsetpcside(kspobject, pc_left, ierr)
        end if
        call echk(ierr, __file__, __line__)
 
        if (trim(kspobjecttype) == 'richardson') then
            call kspsetpcside(kspobject, pc_left, ierr)
            call echk(ierr, __file__, __line__)
        end if
 
        ! Since there is an extra matMult required when using the Richardson preconditioner
        ! with only 1 iteration, only use it when we need to do more than 1 iteration.
        if (globalpreconits > 1) then
            ! Extract preconditioning context for main KSP solver: (master_PC)
            call kspgetpc(kspobject, master_pc, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Set the type of master_PC to ksp. This lets us do multiple
            ! iterations of preconditioner application
            call pcsettype(master_pc, 'ksp', ierr)
            call echk(ierr, __file__, __line__)
 
            ! Get the ksp context from master_PC which is the actual preconditioner:
            call pckspgetksp(master_pc, master_pc_ksp, ierr)
            call echk(ierr, __file__, __line__)
 
            ! master_PC_KSP type will always be of type richardson. If the
            ! number  of iterations is set to 1, this ksp object is transparent.
 
            call kspsettype(master_pc_ksp, 'richardson', ierr)
            call echk(ierr, __file__, __line__)
 
            ! Important to set the norm-type to None for efficiency.
            call kspsetnormtype(master_pc_ksp, ksp_norm_none, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Do one iteration of the outer ksp preconditioners. Note the
            ! tolerances are unsued since we have set KSP_NORM_NON
            call kspsettolerances(master_pc_ksp, petsc_default_real, &
                                  petsc_default_real, petsc_default_real, &
                                  globalpreconits, ierr)
            call echk(ierr, __file__, __line__)
 
            ! Get the 'preconditioner for master_PC_KSP, called 'globalPC'. This
            ! preconditioner is potentially run multiple times.
            call kspgetpc(master_pc_ksp, globalpc, ierr)
            call echk(ierr, __file__, __line__)
        else
            ! Just pull out the pc-object if we are not using kspRichardson
            call kspgetpc(kspobject, globalpc, ierr)
            call echk(ierr, __file__, __line__)
        end if
 
        ! Set the 'globalPC' to additive Schwarz
        call pcsettype(globalpc, 'asm', ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the overlap required
        call pcasmsetoverlap(globalpc, asmoverlap, ierr)
        call echk(ierr, __file__, __line__)
 
        !Setup the main ksp context before extracting the subdomains
        call kspsetup(kspobject, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Extract the ksp objects for each subdomain
        call pcasmgetsubksp(globalpc, nlocal, first, subksp, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Since there is an extra matMult required when using the Richardson preconditioner
        ! with only 1 iteration, only use it when we need to do more than 1 iteration.
        if (localpreconits > 1) then
            ! Set the subksp object to Richardson so we can do multiple iterations on the sub-domains
            call kspsettype(subksp, 'richardson', ierr)
            call echk(ierr, __file__, __line__)
 
            ! Set the number of iterations to do on local blocks. Tolerances are ignored.
            call kspsettolerances(subksp, petsc_default_real, petsc_default_real, petsc_default_real, &
                                  localpreconits, ierr)
            call echk(ierr, __file__, __line__)
 
            ! normtype is NONE because we don't want to check error
            call kspsetnormtype(subksp, ksp_norm_none, ierr)
            call echk(ierr, __file__, __line__)
        else
            ! Set the subksp object to preonly because we are only doing one iteration
            call kspsettype(subksp, 'preonly', ierr)
            call echk(ierr, __file__, __line__)
        end if
 
        ! Extract the preconditioner for subksp object.
        call kspgetpc(subksp, subpc, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the subpc type; only ILU is currently supported
        call pcsettype(subpc, localpctype, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Setup the matrix ordering for the subpc object:
        call pcfactorsetmatorderingtype(subpc, localmatrixordering, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the ILU parameters
        call pcfactorsetlevels(subpc, localfilllevel, ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine setupstandardksp
 
    subroutine setupstandardmultigrid(kspObject, kspObjectType, gmresRestart, preConSide, &
                                      ASMOverlap, outerPreconIts, localMatrixOrdering, fillLevel, localPreConIts, &
                                      ASMOverlapCoarse, fillLevelCoarse, localPreConItsCoarse)
 
        use constants
        use utils, only: echk
        use inputadjoint, only: gmresorthogtype
        use amg, only: amgouterits, amgasmoverlapfine, amgasmoverlapcoarse, amgmatrixordering, &
                       setupshellpc, destroyshellpc, applyshellpc, &
                       amgfilllevelfine, amgfilllevelcoarse, amglocalpreconitsfine, amglocalpreconitscoarse
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Inputs
        ksp kspobject
        character(len=*), intent(in) :: kspObjectType, preConSide, localMatrixOrdering
        integer(kind=intType), intent(in) :: ASMOverlap, fillLevel, gmresRestart, outerPreconIts, localPreConIts
        integer(kind=intType), intent(in) :: ASMOverlapCoarse, fillLevelCoarse, localPreConItsCoarse
 
        ! Working Variables
        pc shellpc
        integer(kind=intType) :: ierr
 
        call kspsettype(kspobject, kspobjecttype, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the preconditioner side from option:
        if (trim(preconside) == 'right') then
            call kspsetpcside(kspobject, pc_right, ierr)
        else
            call kspsetpcside(kspobject, pc_left, ierr)
        end if
        call echk(ierr, __file__, __line__)
 
        call kspgmressetrestart(kspobject, gmresrestart, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Set the orthogonalization method for GMRES
        select case (gmresorthogtype)
        case ('modified_gram_schmidt')
            ! Use modified Gram-Schmidt
            call kspgmressetorthogonalization(kspobject, kspgmresmodifiedgramschmidtorthogonalization, ierr)
        case ('cgs_never_refine')
            ! Use classical Gram-Schmidt with no refinement
            call kspgmressetcgsrefinementtype(kspobject, ksp_gmres_cgs_refine_never, ierr)
        case ('cgs_refine_if_needed')
            ! Use classical Gram-Schmidt with refinement if needed
            call kspgmressetcgsrefinementtype(kspobject, ksp_gmres_cgs_refine_ifneeded, ierr)
        case ('cgs_always_refine')
            ! Use classical Gram-Schmidt with refinement at every iteration
            call kspgmressetcgsrefinementtype(kspobject, ksp_gmres_cgs_refine_always, ierr)
        end select
        call echk(ierr, __file__, __line__)
 
        call kspgetpc(kspobject, shellpc, ierr)
        call echk(ierr, __file__, __line__)
 
        call pcsettype(shellpc, pcshell, ierr)
        call echk(ierr, __file__, __line__)
 
        call pcshellsetsetup(shellpc, setupshellpc, ierr)
        call echk(ierr, __file__, __line__)
 
        call pcshellsetdestroy(shellpc, destroyshellpc, ierr)
        call echk(ierr, __file__, __line__)
 
        call pcshellsetapply(shellpc, applyshellpc, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Save the remaining variables in the AMG module
        amgouterits = outerpreconits
        amgmatrixordering = localmatrixordering
        amgasmoverlapfine = asmoverlap
        amgfilllevelfine = filllevel
        amglocalpreconitsfine = localpreconits
        amgasmoverlapcoarse = asmoverlapcoarse
        amgfilllevelcoarse = filllevelcoarse
        amglocalpreconitscoarse = localpreconitscoarse
 
    end subroutine setupstandardmultigrid
 
    subroutine destroypetscvars
 
        use constants
        use adjointpetsc, only: drdwt, drdwpret, adjointksp, adjointpetscvarsallocated
        use inputadjoint, only: approxpc
        use amg, only: destroyamg
        use utils, only: echk
        implicit none
 
        integer(kind=intType) :: ierr
 
        if (adjointpetscvarsallocated) then
 
            ! Matrices
            call matdestroy(drdwt, ierr)
            call echk(ierr, __file__, __line__)
 
            if (approxpc) then
                call matdestroy(drdwpret, ierr)
                call echk(ierr, __file__, __line__)
            end if
 
            call kspdestroy(adjointksp, ierr)
            call echk(ierr, __file__, __line__)
 
            call destroyamg()
 
            adjointpetscvarsallocated = .false.
        end if
 
    end subroutine destroypetscvars
 
    subroutine initializepetsc
 
        ! Call the C-version of the petsc initialize routine
 
        use adjointpetsc, only: petsc_comm_world
        use communication, only: adflow_comm_world
        implicit none
 
        petsc_comm_world = adflow_comm_world
        call initpetscwrap()
 
    end subroutine initializepetsc
 
    subroutine finalizepetsc
        !
        !      Finalize PETSc by calling the appropriate routine
        !      PetscFinalize provided in the PETSc library. This
        !      automatically calls MPI_Finalize().
        !
        use adjointpetsc, only: petscierr
        implicit none
        call petscfinalize(petscierr)
    end subroutine finalizepetsc
 
    subroutine statepreallocation(onProc, offProc, wSize, stencil, N_stencil, &
                                  level, transposed)
 
        ! This is a generic function that determines the correct
        ! pre-allocation for on and off processor parts of the TRANSPOSED
        ! matrix. With overset, it is quite tricky to determine the
        ! transpose sparsity structure exactly, so we use an alternative
        ! approach. We proceed to determine the non-zeros of the untranposed
        ! matrix, but instead of assigning a non-zero the row we're looping
        ! over, we assign it to the column, which will become a row in the
        ! tranposed matrix. Since this requires communication we use a petsc
        ! vector for doing off processor values. This is not strictly
        ! correct since we will be using the real values as floats, but
        ! since the number of non-zeros per row is always going to be
        ! bounded, we don't have to worry about the integer/floating point
        ! conversions.
 
        use constants
        use blockpointers, only: ndom, il, jl, kl, fringes, flowdoms, globalcell, &
                                 iblank, gind
        use communication, only: adflow_comm_world
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers, echk
        use sorting, only: unique
#include <petsc/finclude/petsc.h>
        use petsc
        implicit none
 
        ! Subroutine Arguments
        integer(kind=intType), intent(in) :: wSize
        integer(kind=intType), intent(in) :: N_stencil
        integer(kind=intType), intent(in) :: stencil(N_stencil, 3)
        integer(kind=intType), intent(out) :: onProc(wSize), offProc(wSize)
        integer(kind=intType), intent(in) :: level
        logical, intent(in) :: transposed
 
        ! Local Variables
        integer(kind=intType) :: nn, i, j, k, sps, ii, jj, kk, iii, jjj, kkk, n, m, gc
        integer(kind=intType) :: iRowStart, iRowEnd, ierr, fInd
        integer(kind=intType), dimension((N_stencil - 1)*8 + 1) :: cellBuffer, dummy
        vec offprocvec
        logical :: overset
        real(kind=realtype), pointer :: tmppointer(:)
 
        call veccreatempi(adflow_comm_world, wsize, petsc_determine, offprocvec, ierr)
        call echk(ierr, __file__, __line__)
 
        ! Zero the cell movement counter
        ii = 0
 
        ! Set the onProc values for each cell to the number of "OFF" time
        ! spectral instances. The "on" spectral instances are accounted for
        ! in the stencil
        onproc(:) = ntimeintervalsspectral - 1
        offproc(:) = 0
        ! Determine the range of onProc in dRdwT
        irowstart = flowdoms(1, 1, 1)%globalCell(2, 2, 2)
        call setpointers(ndom, 1, ntimeintervalsspectral)
        irowend = flowdoms(ndom, 1, ntimeintervalsspectral)%globalCell(il, jl, kl)
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
                ! Loop over each real cell
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
 
                            ! Increment the running ii counter ONLY for each each
                            ! movement of center cell
                            ii = ii + 1
 
                            ! Reset the running tally of the number of neighbours
                            n = 0
 
                            blankedtest: if (iblank(i, j, k) == 1) then
 
                                ! Short-cut flag for cells without interpolated
                                ! cells in it's stencil
                                overset = .false.
 
                                ! Loop over the cells in the provided stencil:
                                do jj = 1, n_stencil
 
                                    ! Determine the cell we are dealing with
                                    iii = stencil(jj, 1) + i
                                    jjj = stencil(jj, 2) + j
                                    kkk = stencil(jj, 3) + k
 
                                    ! Index of the cell we are dealing with. Make
                                    ! code easier to read
                                    gc = globalcell(iii, jjj, kkk)
 
                                    ! Check if the cell in question is a fringe or not:
                                    if (iblank(iii, jjj, kkk) == 1) then
                                        ! regular cell, add to our list, if it is
                                        ! not a boundary
                                        if (gc >= 0) then
                                            n = n + 1
                                            cellbuffer(n) = gc
                                        end if
 
                                    else if (iblank(iii, jjj, kkk) == -1) then
                                        ! Fringe cell. What we do here is loop over
                                        ! the donors for this cell and add any
                                        ! entries that are real cells
                                        overset = .true.
                                        do kk = 1, 8
                                            gc = gind(kk, iii, jjj, kkk)
                                            if (gc >= 0) then
                                                n = n + 1
                                                cellbuffer(n) = gc
                                            end if
                                        end do
                                    end if
                                end do
 
                                ! We have now added 'n' cells to our buffer. For
                                ! the overset interpolation case, it is possible
                                ! (actually highly likely) that the same donor
                                ! cells are used in multiple fringes. To avoid
                                ! allocating more space than necessary, we
                                ! unique-ify the values, producing 'm' unique
                                ! values. If overset wasn't present, we can be
                                ! sure that m=n and we simply don't do the unique
                                ! operation.
 
                                if (overset) then
                                    call unique(cellbuffer, n, m, dummy)
                                else
                                    m = n
                                end if
 
                                ! -------------------- Non-transposed code ----------------
                                if (.not. transposed) then
                                    ! Now we loop over the total number of
                                    ! (unique) neighbours we have and assign them
                                    ! to either an on-proc or an off-proc entry:
                                    do jj = 1, m
                                        gc = cellbuffer(jj)
 
                                        if (gc >= irowstart .and. gc <= irowend) then
                                            onproc(ii) = onproc(ii) + 1
                                        else
                                            offproc(ii) = offproc(ii) + 1
                                        end if
                                    end do
                                else
                                    ! -------------------- Ttransposed code ----------------
 
                                    ! Now we ALSO loop over the total number of
                                    ! (unique) neighbours. However, instead of
                                    ! adding to the non-zeros to the on/offproc for
                                    ! row 'ii', we add them to the column index
                                    ! which will be the row index for the
                                    ! transposed matrix.
                                    do jj = 1, m
                                        gc = cellbuffer(jj)
 
                                        if (gc >= irowstart .and. gc <= irowend) then
                                            ! On processor values can be dealt with
                                            ! directly since the diagonal part is square.
                                            onproc(gc - irowstart + 1) = onproc(gc - irowstart + 1) + 1
                                        else
                                            ! The offproc values need to be sent to
                                            ! the other processors and summed.
                                            call vecsetvalue(offprocvec, gc, real(1), add_values, ierr)
                                            call echk(ierr, __file__, __line__)
                                        end if
                                    end do
                                end if
                            else
                                ! Blanked and interpolated cells only need a single
                                ! non-zero per row for the identity on the diagonal.
                                onproc(ii) = onproc(ii) + 1
                            end if blankedtest
                        end do ! I loop
                    end do ! J loop
                end do ! K loop
            end do ! sps loop
        end do ! Domain Loop
 
        ! Assemble the offproc vector. This doesn't take any time for the
        ! non-transposed operation.
        call vecassemblybegin(offprocvec, ierr)
        call echk(ierr, __file__, __line__)
 
        call vecassemblyend(offprocvec, ierr)
        call echk(ierr, __file__, __line__)
 
        if (transposed) then
            ! Pull the local vector out and convert it back to integers.
            call vecgetarrayf90(offprocvec, tmppointer, ierr)
            call echk(ierr, __file__, __line__)
            do i = 1, wsize
                offproc(i) = int(tmppointer(i) + half) ! Make sure, say 14.99999 is 15.
            end do
 
            call vecrestorearrayf90(offprocvec, tmppointer, ierr)
            call echk(ierr, __file__, __line__)
        end if
 
        ! Done with the temporary offProcVec
        call vecdestroy(offprocvec, ierr)
        call echk(ierr, __file__, __line__)
 
    end subroutine statepreallocation
    subroutine referenceshocksensor
 
        ! Compute the reference shock sensor for PC computations
        use constants
        use blockpointers, only: ib, jb, kb, il, jl, kl, ie, je, ke, shocksensor, &
                                 w, gamma, p, ndom, flowdoms
        use inputtimespectral, only: ntimeintervalsspectral
        use inputphysics, only: equations
        use inputdiscretization, only: spacediscr
        use utils, only: setpointers, echk
        implicit none
 
        ! Working variables
        integer(kind=intType) :: nn, level, sps, i, j, k
 
        level = 1
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
 
                if (equations == eulerequations .or. spacediscr == dissmatrix) then
                    !shockSensor is Pressure
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                shocksensor(i, j, k) = p(i, j, k)
                            end do
                        end do
                    end do
                else
                    ! Enthalpy is used instead
                    do k = 0, kb
                        do j = 2, jl
                            do i = 2, il
                                shocksensor(i, j, k) = p(i, j, k) / (w(i, j, k, irho)**gamma(i, j, k))
                            end do
                        end do
                    end do
 
                    do k = 2, kl
                        do j = 2, jl
                            shocksensor(0, j, k) = p(0, j, k) / (w(0, j, k, irho)**gamma(0, j, k))
                            shocksensor(1, j, k) = p(1, j, k) / (w(1, j, k, irho)**gamma(1, j, k))
                            shocksensor(ie, j, k) = p(ie, j, k) / (w(ie, j, k, irho)**gamma(ie, j, k))
                            shocksensor(ib, j, k) = p(ib, j, k) / (w(ib, j, k, irho)**gamma(ib, j, k))
                        end do
                    end do
 
                    do k = 2, kl
                        do i = 2, il
                            shocksensor(i, 0, k) = p(i, 0, k) / (w(i, 0, k, irho)**gamma(i, 0, k))
                            shocksensor(i, 1, k) = p(i, 1, k) / (w(i, 1, k, irho)**gamma(i, 1, k))
                            shocksensor(i, je, k) = p(i, je, k) / (w(i, je, k, irho)**gamma(i, je, k))
                            shocksensor(i, jb, k) = p(i, jb, k) / (w(i, jb, k, irho)**gamma(i, jb, k))
                        end do
                    end do
                end if
            end do
        end do
    end subroutine referenceshocksensor
 
    subroutine setfdreference(level)
        use constants
        use blockpointers, only: ndom, flowdoms, ib, jb, kb, il, jl, kl, &
                                 shocksensor, w, volref, dw
        use inputphysics, only: liftdirection, veldirfreestream
        use flowvarrefstate, only: nw, nwf
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: echk, setpointers, getdirangle
        use residuals, only: initres_block
        use masterroutines, only: block_res_state
 
        implicit none
 
        ! Input Parameters
        integer(kind=intType) :: level
        ! Working Parameters
        integer(kind=intType) :: i, j, k, l, nn, sps
 
        ! Compute the reference values for doing jacobian with FD
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
 
                call setpointers(nn, level, sps)
                call block_res_state(nn, sps)
                ! Set the values
                do l = 1, nw
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                flowdoms(nn, 1, sps)%wtmp(i, j, k, l) = w(i, j, k, l)
                                flowdoms(nn, 1, sps)%dwtmp(i, j, k, l) = dw(i, j, k, l)
                            end do
                        end do
                    end do
                end do
 
                call initres_block(1, nwf, nn, sps)
 
                ! Note: we have to divide by the volume for dwtmp2 since
                ! normally, dw would have been mulitpiled by 1/Vol in block_res_state
 
                do l = 1, nw
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                flowdoms(nn, 1, sps)%dwtmp2(i, j, k, l) = &
                                    dw(i, j, k, l) / volref(i, j, k)
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end subroutine setfdreference
 
    subroutine resetfdreference(level)
 
        use constants
        use blockpointers, only: ndom, flowdoms, ib, jb, kb, w, dw
        use flowvarrefstate, only: nw, nwf
        use inputtimespectral, only: ntimeintervalsspectral
        use utils, only: setpointers
        implicit none
 
        ! Input Parameters
        integer(kind=intType) :: level
 
        ! Working Parameters
        integer(kind=intType) :: i, j, k, l, nn, sps
        real(kind=realtype) :: sepsensor, cavitation, axismoment
 
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers(nn, level, sps)
                ! Reset w and dw
                do l = 1, nw
                    do k = 0, kb
                        do j = 0, jb
                            do i = 0, ib
                                w(i, j, k, l) = flowdoms(nn, 1, sps)%wtmp(i, j, k, l)
                                dw(i, j, k, l) = flowdoms(nn, 1, sps)%dwtmp(i, j, k, l)
                            end do
                        end do
                    end do
                end do
            end do
        end do
    end subroutine resetfdreference
 
    subroutine setdiffsizes
        !
        !       This routine set the sizes for the pointers that will be
        !       used in the forward debug mode and reverse mode AD.
        !
        use constants
        use blockpointers, only: flowdoms, ib, jb, kb, ie, je, ke, ib, jb, ke, &
                                 nbocos, nviscbocos, ndom
        use flowvarrefstate, only: nw, nwf, nt1, nt2
        use inputtimespectral, only: ntimeintervalsspectral
        use inputphysics, only: equations
        use diffsizes
        implicit none
 
        ! local variables
        integer(kind=intType) :: nLevels
 
        ! Compute nlevels
        nlevels = ubound(flowdoms, 2)
 
        ! Set the size for dynamic pointers to zero for debug purpose
        ! bcdata%norm
        isize1ofdrfdrfbcdata_norm = 0
        isize2ofdrfdrfbcdata_norm = 0
        isize3ofdrfdrfbcdata_norm = 0
 
        ! bcdata%rface
        isize1ofdrfdrfbcdata_rface = 0
        isize2ofdrfdrfbcdata_rface = 0
 
        ! bcdata%m
        isize1ofdrfdrfbcdata_m = 0
        isize2ofdrfdrfbcdata_m = 0
        isize3ofdrfdrfbcdata_m = 0
 
        ! bcdata%fp
        isize1ofdrfdrfbcdata_fp = 0
        isize2ofdrfdrfbcdata_fp = 0
        isize3ofdrfdrfbcdata_fp = 0
 
        ! bcdata%fv
        isize1ofdrfdrfbcdata_fv = 0
        isize2ofdrfdrfbcdata_fv = 0
        isize3ofdrfdrfbcdata_fv = 0
 
        ! bcdata%fp
        isize1ofdrfdrfbcdata_oarea = 0
        isize2ofdrfdrfbcdata_oarea = 0
        isize3ofdrfdrfbcdata_oarea = 0
 
        ! sepSensor
        isize1ofdrfdrfbcdata_sepsensor = 0
        isize2ofdrfdrfbcdata_sepsensor = 0
 
        ! sepSensorKs
        isize1ofdrfdrfbcdata_sepsensorks = 0
        isize2ofdrfdrfbcdata_sepsensorks = 0
 
        ! sepSensorKsArea
        isize1ofdrfdrfbcdata_sepsensorksarea = 0
        isize2ofdrfdrfbcdata_sepsensorksarea = 0
 
        ! Cavitation
        isize1ofdrfdrfbcdata_cavitation = 0
        isize2ofdrfdrfbcdata_cavitation = 0
 
        ! AxisMoment
        isize1ofdrfdrfbcdata_axismoment = 0
        isize2ofdrfdrfbcdata_axismoment = 0
 
        ! viscsubface%tau
        isize1ofdrfdrfviscsubface_tau = 0
        isize2ofdrfdrfviscsubface_tau = 0
        isize3ofdrfdrfviscsubface_tau = 0
 
        ! prod
        isize1ofdrfprod = 0
        isize2ofdrfprod = 0
        isize3ofdrfprod = 0
 
        ! vort
        isize1ofdrfvort = 0
        isize2ofdrfvort = 0
        isize3ofdrfvort = 0
 
        ! dvt
        isize1ofdrfdvt = 0
        isize2ofdrfdvt = 0
        isize3ofdrfdvt = 0
        isize4ofdrfdvt = 0
 
        ! vol
        isize1ofdrfvol = 0
        isize2ofdrfvol = 0
        isize3ofdrfvol = 0
 
        ! rho,etot,u,v,w,p,k
        isize1ofrho = 0
        isize1ofetot = 0
        isize1ofu = 0
        isize1ofv = 0
        isize1ofw = 0
        isize1ofp = 0
        isize1ofk = 0
 
        ! Du1, Du2, Du3
        isize1ofdu1 = 0
        isize1ofdu2 = 0
        isize1ofdu3 = 0
 
        ! Left, Right, Flux
        isize1ofleft = 0
        isize1ofright = 0
        isize1offlux = 0
 
        ! bcdata
        isize1ofdrfbcdata = 0!nbocos
 
        ! s
        isize1ofdrfs = 0!ie
        isize2ofdrfs = je
        isize3ofdrfs = ke
        isize4ofdrfs = 3
 
        ! sfacei
        isize3ofdrfsfacei = 0!ie + 1
        isize2ofdrfsfacei = je
        isize1ofdrfsfacei = ke
 
        ! sfacej
        isize3ofdrfsfacej = 0!ie
        isize2ofdrfsfacej = je + 1
        isize1ofdrfsfacej = ke
 
        ! sfacek
        isize3ofdrfsfacek = 0!ie
        isize2ofdrfsfacek = je
        isize1ofdrfsfacek = ke + 1
 
        ! Define size for the pointers
        ! flowdoms
        isize1ofdrfflowdoms = ndom
        isize2ofdrfflowdoms = nlevels
        isize3ofdrfflowdoms = ntimeintervalsspectral
 
        !viscSubface
        isize1ofdrfviscsubface = nviscbocos
        isize1ofdrfflowdoms_bcdata = nbocos
 
        ! x
        isize4ofdrfx = 3
        isize1ofdrfx = ie + 1
        isize2ofdrfx = je + 1
        isize3ofdrfx = ke + 1
 
        ! flowdoms_x
        isize4ofdrfflowdoms_x = 3
        isize1ofdrfflowdoms_x = ie + 1
        isize2ofdrfflowdoms_x = je + 1
        isize3ofdrfflowdoms_x = ke + 1
 
        if (equations == ransequations) then
            ! rev
            isize1ofdrfrev = ib + 1
            isize2ofdrfrev = jb + 1
            isize3ofdrfrev = kb + 1
        else
            ! rev
            isize1ofdrfrev = 0
            isize2ofdrfrev = 0
            isize3ofdrfrev = 0
        end if
 
        ! rlv
        isize1ofdrfrlv = ib + 1
        isize2ofdrfrlv = jb + 1
        isize3ofdrfrlv = kb + 1
 
        ! w
        isize4ofdrfw = nw
        isize1ofdrfw = ib + 1
        isize2ofdrfw = jb + 1
        isize3ofdrfw = kb + 1
 
        ! flowdoms_x
        isize4ofdrfflowdoms_w = nw
        isize1ofdrfflowdoms_w = ib + 1
        isize2ofdrfflowdoms_w = jb + 1
        isize3ofdrfflowdoms_w = kb + 1
 
        ! flowdoms_dw
        isize4ofdrfflowdoms_dw = nw
        isize1ofdrfflowdoms_dw = ib + 1
        isize2ofdrfflowdoms_dw = jb + 1
        isize3ofdrfflowdoms_dw = kb + 1
 
        ! flowdoms_vol
        isize1ofdrfflowdoms_vol = ib + 1
        isize2ofdrfflowdoms_vol = jb + 1
        isize3ofdrfflowdoms_vol = kb + 1
 
        ! fw
        isize4ofdrffw = nwf
        isize1ofdrffw = ib + 1
        isize2ofdrffw = jb + 1
        isize3ofdrffw = kb + 1
 
        ! dw
        isize4ofdrfdw = nw
        isize1ofdrfdw = ib + 1
        isize2ofdrfdw = jb + 1
        isize3ofdrfdw = kb + 1
 
        ! p
        isize1ofdrfp = ib + 1
        isize2ofdrfp = jb + 1
        isize3ofdrfp = kb + 1
 
        ! gamma
        isize1ofdrfgamma = ib + 1
        isize2ofdrfgamma = jb + 1
        isize3ofdrfgamma = kb + 1
 
        ! dtl
        isize1ofdrfdtl = ie
        isize2ofdrfdtl = je
        isize3ofdrfdtl = ke
 
        ! radI
        isize1ofdrfradi = ie
        isize2ofdrfradi = je
        isize3ofdrfradi = ke
 
        ! radJ
        isize1ofdrfradj = ie
        isize2ofdrfradj = je
        isize3ofdrfradj = ke
 
        ! radK
        isize1ofdrfradk = ie
        isize2ofdrfradk = je
        isize3ofdrfradk = ke
 
        ! sI
        isize1ofdrfsi = ie + 1
        isize2ofdrfsi = je
        isize3ofdrfsi = ke
        isize4ofdrfsi = 3
 
        ! sJ
        isize1ofdrfsj = ie
        isize2ofdrfsj = je + 1
        isize3ofdrfsj = ke
        isize4ofdrfsj = 3
 
        ! sK
        isize1ofdrfsk = ie
        isize2ofdrfsk = je
        isize3ofdrfsk = ke + 1
        isize4ofdrfsk = 3
 
        !bmti1
        isize1ofdrfbmti1 = je
        isize2ofdrfbmti1 = ke
        isize3ofdrfbmti1 = nt2 - nt1 + 1
        isize4ofdrfbmti1 = nt2 - nt1 + 1
 
        !bmti2
        isize1ofdrfbmti2 = je
        isize2ofdrfbmti2 = ke
        isize3ofdrfbmti2 = nt2 - nt1 + 1
        isize4ofdrfbmti2 = nt2 - nt1 + 1
 
        !bmtj1
        isize1ofdrfbmtj1 = ie
        isize2ofdrfbmtj1 = ke
        isize3ofdrfbmtj1 = nt2 - nt1 + 1
        isize4ofdrfbmtj1 = nt2 - nt1 + 1
 
        !bmtj2
        isize1ofdrfbmtj2 = ie
        isize2ofdrfbmtj2 = ke
        isize3ofdrfbmtj2 = nt2 - nt1 + 1
        isize4ofdrfbmtj2 = nt2 - nt1 + 1
 
        !bmtk1
        isize1ofdrfbmtk1 = ie
        isize2ofdrfbmtk1 = je
        isize3ofdrfbmtk1 = nt2 - nt1 + 1
        isize4ofdrfbmtk1 = nt2 - nt1 + 1
 
        !bmtk2
        isize1ofdrfbmtk2 = ie
        isize2ofdrfbmtk2 = je
        isize3ofdrfbmtk2 = nt2 - nt1 + 1
        isize4ofdrfbmtk2 = nt2 - nt1 + 1
 
        !bvti1
        isize1ofdrfbvti1 = je
        isize2ofdrfbvti1 = ke
        isize3ofdrfbvti1 = nt2 - nt1 + 1
 
        !bvti2
        isize1ofdrfbvti2 = je
        isize2ofdrfbvti2 = ke
        isize3ofdrfbvti2 = nt2 - nt1 + 1
        !bvti1
        isize1ofdrfbvti1 = je
        isize2ofdrfbvti1 = ke
        isize3ofdrfbvti1 = nt2 - nt1 + 1
 
        !bvti2
        isize1ofdrfbvti2 = je
        isize2ofdrfbvti2 = ke
        isize3ofdrfbvti2 = nt2 - nt1 + 1
 
        !bvtk1
        isize1ofdrfbvtk1 = ie
        isize2ofdrfbvtk1 = je
        isize3ofdrfbvtk1 = nt2 - nt1 + 1
 
        !bvtk2
        isize1ofdrfbvtk2 = ie
        isize2ofdrfbvtk2 = je
        isize3ofdrfbvtk2 = nt2 - nt1 + 1
 
    end subroutine setdiffsizes
 
end module adjointutils