adjointDebug.F90

Go to the documentation of this file.

! This module is used for debugging and testing only
 
module adjointdebug
 
contains
 
#ifndef USE_COMPLEX
 
    subroutine computematrixfreeproductfwdfd(xvdot, extradot, wdot, bcDataValuesdot, &
                                             useSpatial, useState, famLists, &
                                             bcDataNames, bcDataValues, bcDataFamLists, bcVarsEmpty, &
                                             dwdot, funcsDot, fDot, &
                                             costSize, fSize, nTime, h)
 
        ! This routine is used to debug master_d. It uses the forward seeds to set perturbations
        ! and then computes the value of the derivatives using forward finite diffenece
 
        use constants
        use adjointvars
        use blockpointers, only: ndom
        use communication, only: adflow_comm_world
        use inputtimespectral, only: ntimeintervalsspectral
        use inputphysics, only: pointref, alpha, beta, equations, machcoef, &
                                mach, machgrid, rgasdim
        use iteration, only: currentlevel, groundlevel
        use flowvarrefstate, only: pinfdim, rhoinfdim, tinfdim
        use blockpointers, only: ndom, il, jl, kl, wd, x, w, dw, dwd, nbocos, nviscbocos
 
        use adjointutils, only: allocderivativevalues, zeroadseeds
        use masterroutines, only: master
        use utils, only: iswalltype, setpointers, setpointers_d, echk
        use flowvarrefstate, only: nw, nwf
        use walldistancedata, only: xsurf, xsurfvec
        use walldistance, only: updatexsurf
        implicit none
        !
        ! Input Variables
        !
        ! derivative seeds
        real(kind=realtype), dimension(:), intent(in) :: xvdot
        real(kind=realtype), dimension(:), intent(in) :: extradot
        real(kind=realtype), dimension(:), intent(in) :: wdot
 
        ! size data
        logical, intent(in) :: useSpatial, useState
        integer(kind=intType), dimension(:, :) :: famLists
        integer(kind=intType) :: costSize, fSize, nTime
 
        character, dimension(:, :), intent(in) :: bcDataNames
        real(kind=realtype), dimension(:), intent(inout) :: bcdatavalues, bcdatavaluesdot
        integer(kind=intType), dimension(:, :) :: bcDataFamLists
        logical, intent(in) :: BCVarsEmpty
 
        ! Finite difference parameters
        real(kind=realtype), intent(in) :: h ! step size for Finite Difference
 
        !
        ! Ouput Variables
        !
        ! Output derivative seeds
        real(kind=realtype), dimension(size(wDot)), intent(out) :: dwdot
        real(kind=realtype), dimension(costSize, size(famLists, 1)), intent(out) :: funcsdot
        real(kind=realtype), dimension(3, fSize, nTime), intent(out) :: fdot
 
        !
        ! Working Variables
        !
        integer(kind=intType) :: nn, sps, level
        integer(kind=intType) :: ierr, mm, i, j, k, l, ii, jj, iRegion
 
        real(kind=realtype), dimension(costSize, size(famLists, 1)) :: funcs
 
        ! Input Arguments for master:
        real(kind=realtype), dimension(costSize, size(famLists, 1)) :: funcvalues
        real(kind=realtype), dimension(:, :, :), allocatable :: forces
 
        fsize = size(fdot, 2)
        allocate (forces(3, fsize, ntimeintervalsspectral))
 
        ! Need to trick the residual evalution to use coupled (mean flow and
        ! turbulent) together.
        level = 1
        currentlevel = level
        groundlevel = level
 
        ! Allocate the memory we need for derivatives if not done so
        ! already. Note this isn't deallocated until the adflow is
        ! destroyed.
 
        ! a bit over kill for our needs, but so values are stored in the seeds
        if (.not. derivvarsallocated) then
            call allocderivativevalues(level)
        end if
 
        ! Zero all AD seesd.
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call zeroadseeds(nn, level, sps)
            end do
        end do
 
        ! ----------------------------- Run Master ---------------------------------
        ! Run the super-dee-duper master rotuine
        if (bcvarsempty) then
            call master(usespatial, &
                        famlists, funcvalues, &
                        forces)
 
        else
            call master(usespatial, &
                        famlists, funcvalues, &
                        forces, &
                        bcdatanames, bcdatavalues, bcdatafamlists)
        end if
 
        ! Copy base val (f) into variables for the final vals (f(x+dx) - f)/dx
        ! we add the negative sign here instead of doing it later
 
        ii = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                ii = ii + 1
                                dwd(i, j, k, l) = -dw(i, j, k, l)
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
        fdot = -forces
        funcsdot = -funcvalues
 
        !  --------------------- apply the perturbations ----------------------------
        ! Set the extra seeds now do the extra ones. Note that we are assuming the
        ! machNumber used for the coefficients follows the Mach number,
        ! not the grid mach number.
 
        alpha = alpha + h * extradot(ialpha)
        beta = beta + h * extradot(ibeta)
        mach = mach + h * extradot(imach)
        machcoef = machcoef + h * extradot(imach)
        machgrid = machgrid + h * extradot(imachgrid)
        pinfdim = pinfdim + h * extradot(ipressure)
        rhoinfdim = rhoinfdim + h * extradot(idensity)
        tinfdim = tinfdim + h * extradot(itemperature)
        pointref(1) = pointref(1) + h * extradot(ipointrefx)
        pointref(2) = pointref(2) + h * extradot(ipointrefy)
        pointref(3) = pointref(3) + h * extradot(ipointrefz)
        rgasdim = rgasdim + h * zero
 
        ! Set the provided w and x seeds:
        ii = 0
        jj = 0
        domainloop1: do nn = 1, ndom
            spectalloop1: do sps = 1, ntimeintervalsspectral
                call setpointers(nn, 1, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                ii = ii + 1
                                x(i, j, k, l) = x(i, j, k, l) + xvdot(ii) * h
                            end do
                        end do
                    end do
                end do
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                jj = jj + 1
                                w(i, j, k, l) = w(i, j, k, l) + wdot(jj) * h
                            end do
                        end do
                    end do
                end do
            end do spectalloop1
        end do domainloop1
 
        bcdatavalues = bcdatavalues + bcdatavaluesdot * h
 
        ! The xvolume that is set by used by update Xsurf is only allocated with
        ! rans
        if (equations == ransequations) then
            call updatexsurf(level)
        end if
 
        ! ----------------------------- Run Master ---------------------------------
        ! Run the super-dee-duper master rotuine
        if (bcvarsempty) then
            call master(usespatial, &
                        famlists, funcvalues, &
                        forces)
 
        else
            call master(usespatial, &
                        famlists, funcvalues, &
                        forces, &
                        bcdatanames, bcdatavalues, bcdatafamlists)
 
        end if
 
        ! ------------------------- set the output vectors ----------------------
        ! Set the extra seeds now do the extra ones. Note that we are assuming the
        ! machNumber used for the coefficients follows the Mach number,
        ! not the grid mach number.
 
        alpha = alpha - h * extradot(ialpha)
        beta = beta - h * extradot(ibeta)
        mach = mach - h * extradot(imach)
        machcoef = machcoef - h * extradot(imach)
        machgrid = machgrid - h * extradot(imachgrid)
        pinfdim = pinfdim - h * extradot(ipressure)
        rhoinfdim = rhoinfdim - h * extradot(idensity)
        tinfdim = tinfdim - h * extradot(itemperature)
        pointref(1) = pointref(1) - h * extradot(ipointrefx)
        pointref(2) = pointref(2) - h * extradot(ipointrefy)
        pointref(3) = pointref(3) - h * extradot(ipointrefz)
        rgasdim = rgasdim - h * zero
 
        ! Copy out the residual derivative into the provided dwDot and remove the
        ! perturbation
        ii = 0
        jj = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                ii = ii + 1
                                x(i, j, k, l) = x(i, j, k, l) - xvdot(ii) * h
                            end do
                        end do
                    end do
                end do
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                jj = jj + 1
                                w(i, j, k, l) = w(i, j, k, l) - wdot(jj) * h
                                dwd(i, j, k, l) = (dwd(i, j, k, l) + dw(i, j, k, l)) / h
                                dwdot(jj) = dwd(i, j, k, l) ! copy values to output
                            end do
                        end do
                    end do
                end do
 
            end do
        end do
 
        if (equations == ransequations) then
            call updatexsurf(level)
        end if
        bcdatavalues = bcdatavalues - bcdatavaluesdot * h
 
        fdot = (fdot + forces) / h
        funcsdot = (funcsdot + funcvalues) / h
 
    end subroutine computematrixfreeproductfwdfd
 
    subroutine printadseeds(nn, level, sps)
!DIR$ NOOPTIMIZE
        ! this routine is used for debugging master_d, and master_b.
        ! it prints all the AD seeds used
        ! it is useful to save the output and compare it with a diff tool
 
        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
        write (*, *) 'ADSeeds for block', nn, ' at level ', level, ' at sps ', sps
        write (*, *) 'd2wall ', minval(flowdomsd(nn, level, sps)%d2wall), &
            maxval(flowdomsd(nn, level, sps)%d2wall), &
            norm2(flowdomsd(nn, level, sps)%d2wall)
        write (*, *) 'x ', minval(flowdomsd(nn, level, sps)%x), &
            maxval(flowdomsd(nn, level, sps)%x), &
            norm2(flowdomsd(nn, level, sps)%x)
        write (*, *) 'si ', minval(flowdomsd(nn, level, sps)%si), &
            maxval(flowdomsd(nn, level, sps)%si), &
            norm2(flowdomsd(nn, level, sps)%si)
        write (*, *) 'sj ', minval(flowdomsd(nn, level, sps)%sj), &
            maxval(flowdomsd(nn, level, sps)%sj), &
            norm2(flowdomsd(nn, level, sps)%sj)
        write (*, *) 'sk ', minval(flowdomsd(nn, level, sps)%sk), &
            maxval(flowdomsd(nn, level, sps)%sk), &
            norm2(flowdomsd(nn, level, sps)%sk)
        write (*, *) 'vol ', minval(flowdomsd(nn, level, sps)%vol), &
            maxval(flowdomsd(nn, level, sps)%vol), &
            norm2(flowdomsd(nn, level, sps)%vol)
 
        write (*, *) 's ', minval(flowdomsd(nn, level, sps)%s), &
            maxval(flowdomsd(nn, level, sps)%s), &
            norm2(flowdomsd(nn, level, sps)%s)
        write (*, *) 'sFaceI ', minval(flowdomsd(nn, level, sps)%sFaceI), &
            maxval(flowdomsd(nn, level, sps)%sFaceI), &
            norm2(flowdomsd(nn, level, sps)%sFaceI)
        write (*, *) 'sFaceJ ', minval(flowdomsd(nn, level, sps)%sFaceJ), &
            maxval(flowdomsd(nn, level, sps)%sFaceJ), &
            norm2(flowdomsd(nn, level, sps)%sFaceJ)
        write (*, *) 'sFaceK ', minval(flowdomsd(nn, level, sps)%sFaceK), &
            maxval(flowdomsd(nn, level, sps)%sFaceK), &
            norm2(flowdomsd(nn, level, sps)%sFaceK)
 
        write (*, *) 'w ', minval(flowdomsd(nn, level, sps)%w), &
            maxval(flowdomsd(nn, level, sps)%w), &
            norm2(flowdomsd(nn, level, sps)%w)
        write (*, *) 'dw ', minval(flowdomsd(nn, level, sps)%dw), &
            maxval(flowdomsd(nn, level, sps)%dw), &
            norm2(flowdomsd(nn, level, sps)%dw)
        write (*, *) 'fw ', minval(flowdomsd(nn, level, sps)%fw), &
            maxval(flowdomsd(nn, level, sps)%fw), &
            norm2(flowdomsd(nn, level, sps)%fw)
        write (*, *) 'scratch ', minval(flowdomsd(nn, level, sps)%scratch), &
            maxval(flowdomsd(nn, level, sps)%scratch), &
            norm2(flowdomsd(nn, level, sps)%scratch)
 
        write (*, *) 'p ', minval(flowdomsd(nn, level, sps)%p), &
            maxval(flowdomsd(nn, level, sps)%p), &
            norm2(flowdomsd(nn, level, sps)%p)
        write (*, *) 'gamma ', minval(flowdomsd(nn, level, sps)%gamma), &
            maxval(flowdomsd(nn, level, sps)%gamma), &
            norm2(flowdomsd(nn, level, sps)%gamma)
        write (*, *) 'aa ', minval(flowdomsd(nn, level, sps)%aa), &
            maxval(flowdomsd(nn, level, sps)%aa), &
            norm2(flowdomsd(nn, level, sps)%aa)
 
        write (*, *) 'rlv ', minval(flowdomsd(nn, level, sps)%rlv), &
            maxval(flowdomsd(nn, level, sps)%rlv), &
            norm2(flowdomsd(nn, level, sps)%rlv)
        write (*, *) 'rev ', minval(flowdomsd(nn, level, sps)%rev), &
            maxval(flowdomsd(nn, level, sps)%rev), &
            norm2(flowdomsd(nn, level, sps)%rev)
 
        write (*, *) 'radI ', minval(flowdomsd(nn, level, sps)%radI), &
            maxval(flowdomsd(nn, level, sps)%radI), &
            norm2(flowdomsd(nn, level, sps)%radI)
        write (*, *) 'radJ ', minval(flowdomsd(nn, level, sps)%radJ), &
            maxval(flowdomsd(nn, level, sps)%radJ), &
            norm2(flowdomsd(nn, level, sps)%radJ)
        write (*, *) 'radK ', minval(flowdomsd(nn, level, sps)%radK), &
            maxval(flowdomsd(nn, level, sps)%radK), &
            norm2(flowdomsd(nn, level, sps)%radK)
 
        write (*, *) 'ux ', minval(flowdomsd(nn, level, sps)%ux), &
            maxval(flowdomsd(nn, level, sps)%ux), &
            norm2(flowdomsd(nn, level, sps)%ux)
        write (*, *) 'uy ', minval(flowdomsd(nn, level, sps)%uy), &
            maxval(flowdomsd(nn, level, sps)%uy), &
            norm2(flowdomsd(nn, level, sps)%uy)
        write (*, *) 'uz ', minval(flowdomsd(nn, level, sps)%uz), &
            maxval(flowdomsd(nn, level, sps)%uz), &
            norm2(flowdomsd(nn, level, sps)%uz)
        write (*, *) 'vx ', minval(flowdomsd(nn, level, sps)%vx), &
            maxval(flowdomsd(nn, level, sps)%vx), &
            norm2(flowdomsd(nn, level, sps)%vx)
        write (*, *) 'vy ', minval(flowdomsd(nn, level, sps)%vy), &
            maxval(flowdomsd(nn, level, sps)%vy), &
            norm2(flowdomsd(nn, level, sps)%vy)
        write (*, *) 'vz ', minval(flowdomsd(nn, level, sps)%vz), &
            maxval(flowdomsd(nn, level, sps)%vz), &
            norm2(flowdomsd(nn, level, sps)%vz)
        write (*, *) 'wx ', minval(flowdomsd(nn, level, sps)%wx), &
            maxval(flowdomsd(nn, level, sps)%wx), &
            norm2(flowdomsd(nn, level, sps)%wx)
        write (*, *) 'wy ', minval(flowdomsd(nn, level, sps)%wy), &
            maxval(flowdomsd(nn, level, sps)%wy), &
            norm2(flowdomsd(nn, level, sps)%wy)
        write (*, *) 'wz ', minval(flowdomsd(nn, level, sps)%wz), &
            maxval(flowdomsd(nn, level, sps)%wz), &
            norm2(flowdomsd(nn, level, sps)%wz)
        write (*, *) 'qx ', minval(flowdomsd(nn, level, sps)%qx), &
            maxval(flowdomsd(nn, level, sps)%qx), &
            norm2(flowdomsd(nn, level, sps)%qx)
        write (*, *) 'qy ', minval(flowdomsd(nn, level, sps)%qy), &
            maxval(flowdomsd(nn, level, sps)%qy), &
            norm2(flowdomsd(nn, level, sps)%qy)
        write (*, *) 'qz ', minval(flowdomsd(nn, level, sps)%qz), &
            maxval(flowdomsd(nn, level, sps)%qz), &
            norm2(flowdomsd(nn, level, sps)%qz)
 
        write (*, *) 'bmti1 ', minval(flowdomsd(nn, level, sps)%bmti1), &
            maxval(flowdomsd(nn, level, sps)%bmti1), &
            norm2(flowdomsd(nn, level, sps)%bmti1)
        write (*, *) 'bmti2 ', minval(flowdomsd(nn, level, sps)%bmti2), &
            maxval(flowdomsd(nn, level, sps)%bmti2), &
            norm2(flowdomsd(nn, level, sps)%bmti2)
        write (*, *) 'bmtj1 ', minval(flowdomsd(nn, level, sps)%bmtj1), &
            maxval(flowdomsd(nn, level, sps)%bmtj1), &
            norm2(flowdomsd(nn, level, sps)%bmtj1)
        write (*, *) 'bmtj2 ', minval(flowdomsd(nn, level, sps)%bmtj2), &
            maxval(flowdomsd(nn, level, sps)%bmtj2), &
            norm2(flowdomsd(nn, level, sps)%bmtj2)
        write (*, *) 'bmtk1 ', minval(flowdomsd(nn, level, sps)%bmtk1), &
            maxval(flowdomsd(nn, level, sps)%bmtk1), &
            norm2(flowdomsd(nn, level, sps)%bmtk1)
        write (*, *) 'bmtk2 ', minval(flowdomsd(nn, level, sps)%bmtk2), &
            maxval(flowdomsd(nn, level, sps)%bmtk2), &
            norm2(flowdomsd(nn, level, sps)%bmtk2)
        write (*, *) 'bvti1 ', minval(flowdomsd(nn, level, sps)%bvti1), &
            maxval(flowdomsd(nn, level, sps)%bvti1), &
            norm2(flowdomsd(nn, level, sps)%bvti1)
        write (*, *) 'bvti2 ', minval(flowdomsd(nn, level, sps)%bvti2), &
            maxval(flowdomsd(nn, level, sps)%bvti2), &
            norm2(flowdomsd(nn, level, sps)%bvti2)
        write (*, *) 'bvtj1 ', minval(flowdomsd(nn, level, sps)%bvtj1), &
            maxval(flowdomsd(nn, level, sps)%bvtj1), &
            norm2(flowdomsd(nn, level, sps)%bvtj1)
        write (*, *) 'bvtj2 ', minval(flowdomsd(nn, level, sps)%bvtj2), &
            maxval(flowdomsd(nn, level, sps)%bvtj2), &
            norm2(flowdomsd(nn, level, sps)%bvtj2)
        write (*, *) 'bvtk1 ', minval(flowdomsd(nn, level, sps)%bvtk1), &
            maxval(flowdomsd(nn, level, sps)%bvtk1), &
            norm2(flowdomsd(nn, level, sps)%bvtk1)
        write (*, *) 'bvtk2 ', minval(flowdomsd(nn, level, sps)%bvtk2), &
            maxval(flowdomsd(nn, level, sps)%bvtk2), &
            norm2(flowdomsd(nn, level, sps)%bvtk2)
 
        bocoloop: do mm = 1, flowdoms(nn, level, sps)%nBocos
            write (*, *) 'mm', mm, 'BCData(mm)%norm', minval(flowdomsd(nn, level, sps)%BCData(mm)%norm), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%norm), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%norm)
            write (*, *) 'mm', mm, 'bcData(mm)%rface ', minval(flowdomsd(nn, level, sps)%bcData(mm)%rface), &
                maxval(flowdomsd(nn, level, sps)%bcData(mm)%rface), &
                norm2(flowdomsd(nn, level, sps)%bcData(mm)%rface)
            write (*, *) 'mm', mm, 'bcData(mm)%Fv ', minval(flowdomsd(nn, level, sps)%bcData(mm)%Fv), &
                maxval(flowdomsd(nn, level, sps)%bcData(mm)%Fv), &
                norm2(flowdomsd(nn, level, sps)%bcData(mm)%Fv)
            write (*, *) 'mm', mm, 'bcData(mm)%Fp ', minval(flowdomsd(nn, level, sps)%bcData(mm)%Fp), &
                maxval(flowdomsd(nn, level, sps)%bcData(mm)%Fp), &
                norm2(flowdomsd(nn, level, sps)%bcData(mm)%Fp)
            write (*, *) 'mm', mm, 'bcData(mm)%Tv ', minval(flowdomsd(nn, level, sps)%bcData(mm)%Tv), &
                maxval(flowdomsd(nn, level, sps)%bcData(mm)%Tv), &
                norm2(flowdomsd(nn, level, sps)%bcData(mm)%Tv)
            write (*, *) 'mm', mm, 'bcData(mm)%Tp ', minval(flowdomsd(nn, level, sps)%bcData(mm)%Tp), &
                maxval(flowdomsd(nn, level, sps)%bcData(mm)%Tp), &
                norm2(flowdomsd(nn, level, sps)%bcData(mm)%Tp)
            write (*, *) 'mm', mm, 'bcData(mm)%area ', minval(flowdomsd(nn, level, sps)%bcData(mm)%area), &
                maxval(flowdomsd(nn, level, sps)%bcData(mm)%area), &
                norm2(flowdomsd(nn, level, sps)%bcData(mm)%area)
            write (*, *) 'mm', mm, 'BCData(mm)%uSlip ', minval(flowdomsd(nn, level, sps)%BCData(mm)%uSlip), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%uSlip), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%uSlip)
            write (*, *) 'mm', mm, 'BCData(mm)%TNS_Wall ', minval(flowdomsd(nn, level, sps)%BCData(mm)%TNS_Wall), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%TNS_Wall), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%TNS_Wall)
            write (*, *) 'mm', mm, 'BCData(mm)%ptInlet ', minval(flowdomsd(nn, level, sps)%BCData(mm)%ptInlet), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%ptInlet), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%ptInlet)
            write (*, *) 'mm', mm, 'BCData(mm)%htInlet ', minval(flowdomsd(nn, level, sps)%BCData(mm)%htInlet), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%htInlet), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%htInlet)
            write (*, *) 'mm', mm, 'BCData(mm)%ttInlet ', minval(flowdomsd(nn, level, sps)%BCData(mm)%ttInlet), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%ttInlet), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%ttInlet)
            write (*, *) 'mm', mm, 'BCData(mm)%turbInlet ', minval(flowdomsd(nn, level, sps)%BCData(mm)%turbInlet), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%turbInlet), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%turbInlet)
            write (*, *) 'mm', mm, 'BCData(mm)%ps ', minval(flowdomsd(nn, level, sps)%BCData(mm)%ps), &
                maxval(flowdomsd(nn, level, sps)%BCData(mm)%ps), &
                norm2(flowdomsd(nn, level, sps)%BCData(mm)%ps)
 
        end do bocoloop
 
        viscbocoloop: do mm = 1, flowdoms(nn, level, sps)%nViscBocos
            write (*, *) 'mm', mm, 'viscSubface(mm)%tau ', minval(flowdomsd(nn, level, sps)%viscSubface(mm)%tau), &
                maxval(flowdomsd(nn, level, sps)%viscSubface(mm)%tau), &
                norm2(flowdomsd(nn, level, sps)%viscSubface(mm)%tau)
            write (*, *) 'mm', mm, 'viscSubface(mm)%q ', minval(flowdomsd(nn, level, sps)%viscSubface(mm)%q), &
                maxval(flowdomsd(nn, level, sps)%viscSubface(mm)%q), &
                norm2(flowdomsd(nn, level, sps)%viscSubface(mm)%q)
        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
                write (*, *) 'commPatternOverset(level, sps)%sendList(i)%interpd ', &
                    minval(commpatternoverset(level, sps)%sendList(i)%interpd), &
                    maxval(commpatternoverset(level, sps)%sendList(i)%interpd), &
                    norm2(commpatternoverset(level, sps)%sendList(i)%interpd)
            end do sends
            write (*, *) 'internalOverset(level, sps)%donorInterpd ', &
                minval(internaloverset(level, sps)%donorInterpd), &
                maxval(internaloverset(level, sps)%donorInterpd), &
                norm2(internaloverset(level, sps)%donorInterpd)
        end if
 
        write (*, *) 'alphad ', alphad
        write (*, *) 'betad ', betad
        write (*, *) 'machd ', machd
        write (*, *) 'machGridd ', machgridd
        write (*, *) 'machCoefd ', machcoefd
        write (*, *) 'pinfdimd ', pinfdimd
        write (*, *) 'tinfdimd ', tinfdimd
        write (*, *) 'rhoinfdimd ', rhoinfdimd
        write (*, *) 'rgasdimd ', rgasdimd
        write (*, *) 'pointrefd ', pointrefd
        write (*, *) 'prefd ', prefd
        write (*, *) 'rhoRefd ', rhorefd
        write (*, *) 'Trefd ', trefd
        write (*, *) 'murefd ', murefd
        write (*, *) 'urefd ', urefd
        write (*, *) 'hrefd ', hrefd
        write (*, *) 'timerefd ', timerefd
        write (*, *) 'pinfd ', pinfd
        write (*, *) 'pinfCorrd ', pinfcorrd
        write (*, *) 'rhoinfd ', rhoinfd
        write (*, *) 'uinfd ', uinfd
        write (*, *) 'rgasd ', rgasd
        write (*, *) 'muinfd ', muinfd
        write (*, *) 'gammainfd ', gammainfd
        write (*, *) 'winfd ', winfd
        write (*, *) 'veldirfreestreamd ', veldirfreestreamd
        write (*, *) 'liftdirectiond ', liftdirectiond
        write (*, *) 'dragdirectiond ', dragdirectiond
 
        ! 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)
                                write (*, *) idom, iboco, idata, idirichlet, 'dataArr(:) ' &
                                    , cgnsdomsd(idom)%bocoInfo(iboco)%dataSet(idata)% &
                                    dirichletarrays(idirichlet)%dataArr(:)
                            end do
                        end if
                    end do
                end if
            end do
        end do
 
        ! And the reverse seeds in the actuator zones
        do i = 1, nactuatorregions
            write (*, *) 'actuatorRegionsd(i)%Force ', actuatorregionsd(i)%force
            write (*, *) 'actuatorRegionsd(i)%Torque ', actuatorregionsd(i)%torque
        end do
 
    end subroutine printadseeds
 
#else
 
    subroutine computematrixfreeproductfwdcs(xvdot, extradot, wdot, bcDataValuesdot, &
                                             useSpatial, useState, famLists, &
                                             bcDataNames, bcDataValues, bcDataFamLists, bcVarsEmpty, &
                                             dwdot, funcsDot, fDot, &
                                             costSize, fSize, nTime, h_mag)
 
        ! This routine is used to debug master_d. It uses the forward seeds to set perturbations
        ! and then computes the value of the derivatives using forward complex step
 
        use constants
        use adjointvars
        use blockpointers, only: ndom
        use communication, only: adflow_comm_world
        use inputtimespectral, only: ntimeintervalsspectral
        use inputphysics, only: pointref, alpha, beta, equations, machcoef, &
                                mach, machgrid, rgasdim
        use iteration, only: currentlevel, groundlevel
        use flowvarrefstate, only: pinfdim, rhoinfdim, tinfdim
        use blockpointers, only: ndom, il, jl, kl, wd, x, w, dw, dwd, nbocos, nviscbocos
 
        use adjointutils, only: allocderivativevalues, zeroadseeds
        use masterroutines, only: master
        use utils, only: iswalltype, setpointers, setpointers_d, echk
        use flowvarrefstate, only: nw, nwf
        use walldistancedata, only: xsurf, xsurfvec
        use walldistance, only: updatexsurf
 
        implicit none
 
        ! Input Variables
        complex(kind=realType), dimension(:), intent(in) :: xvdot
        complex(kind=realType), dimension(:), intent(in) :: extradot
        complex(kind=realType), dimension(:), intent(in) :: wdot
 
        logical, intent(in) :: useSpatial, useState
        integer(kind=intType), dimension(:, :) :: famLists
        integer(kind=intType) :: costSize, fSize, nTime
 
        character, dimension(:, :), intent(in) :: bcDataNames
        real(kind=realtype), dimension(:), intent(inout) :: bcdatavalues, bcdatavaluesdot
        integer(kind=intType), dimension(:, :) :: bcDataFamLists
        logical, intent(in) :: BCVarsEmpty
 
        ! step parameters
        real(kind=alwaysrealtype), intent(in) :: h_mag ! step size for step
 
        ! Ouput Variables
        complex(kind=realType), dimension(size(wdot)), intent(out) :: dwDot
        complex(kind=realType), dimension(costSize, size(famLists, 1)), intent(out) :: funcsDot
        complex(kind=realType), dimension(3, fSize, nTime), intent(out) :: fDot
 
        ! Working Variables
        integer(kind=intType) :: nn, sps, level
        integer(kind=intType) :: ierr, mm, i, j, k, l, ii, jj, iRegion
 
        complex(kind=realType), dimension(costSize, size(famLists, 1)) :: funcs
 
        ! Input Arguments for master:
        complex(kind=realType), dimension(costSize, size(famLists, 1)) :: funcValues
 
        ! Working Variables
        complex(kind=realType), dimension(:, :, :), allocatable :: forces
        complex(kind=realType) :: h ! step size for Finite Difference
 
        ! note that h_mag does not have to be complex
        ! it is just the magnitude of the complex perturbation
        h = cmplx(0, h_mag)
 
        fsize = size(fdot, 2)
        allocate (forces(3, fsize, ntimeintervalsspectral))
 
        ! Need to trick the residual evalution to use coupled (mean flow and
        ! turbulent) together.
        level = 1
        currentlevel = level
        groundlevel = level
 
        ! Allocate the memory we need for derivatives if not done so
        ! already. Note this isn't deallocated until the adflow is
        ! destroyed.
        if (.not. derivvarsallocated) then
            call allocderivativevalues(level)
        end if
 
        ! Zero all AD seesd.
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call zeroadseeds(nn, level, sps)
            end do
        end do
 
        ! Set the extra seeds now do the extra ones. Note that we are assuming the
        ! machNumber used for the coefficients follows the Mach number,
        ! not the grid mach number.
 
        alpha = alpha + h * extradot(ialpha)
        beta = beta + h * extradot(ibeta)
        mach = mach + h * extradot(imach)
        machcoef = machcoef + h * extradot(imach)
        machgrid = machgrid + h * extradot(imachgrid)
        pinfdim = pinfdim + h * extradot(ipressure)
        rhoinfdim = rhoinfdim + h * extradot(idensity)
        tinfdim = tinfdim + h * extradot(itemperature)
        pointref(1) = pointref(1) + h * extradot(ipointrefx)
        pointref(2) = pointref(2) + h * extradot(ipointrefy)
        pointref(3) = pointref(3) + h * extradot(ipointrefz)
        rgasdim = rgasdim + h * zero
 
        !  --------------------- apply the perturbations ----------------------------
        ! Set the provided w and x seeds:
        ii = 0
        jj = 0
        domainloop1: do nn = 1, ndom
            spectalloop1: do sps = 1, ntimeintervalsspectral
                call setpointers(nn, 1, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                ii = ii + 1
                                x(i, j, k, l) = x(i, j, k, l) + xvdot(ii) * h
                            end do
                        end do
                    end do
                end do
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                jj = jj + 1
                                w(i, j, k, l) = w(i, j, k, l) + wdot(jj) * h
                            end do
                        end do
                    end do
                end do
            end do spectalloop1
        end do domainloop1
 
        bcdatavalues = bcdatavalues + bcdatavaluesdot * h
 
        if (equations == ransequations) then
            call updatexsurf(level)
        end if
 
        ! ----------------------------- Run Master ---------------------------------
        ! Run the super-dee-duper master rotuine
        if (bcvarsempty) then
            call master(usespatial, &
                        famlists, funcvalues, &
                        forces)
        else
            call master(usespatial, &
                        famlists, funcvalues, &
                        forces, &
                        bcdatanames, bcdatavalues, bcdatafamlists)
        end if
 
        ! Copy out the residual derivative into the provided dwDot and remove the
        ! perturbation
        ii = 0
        jj = 0
        do nn = 1, ndom
            do sps = 1, ntimeintervalsspectral
                call setpointers_d(nn, 1, sps)
                do k = 1, kl
                    do j = 1, jl
                        do i = 1, il
                            do l = 1, 3
                                ii = ii + 1
                                x(i, j, k, l) = x(i, j, k, l) - xvdot(ii) * h
                            end do
                        end do
                    end do
                end do
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            do l = 1, nw
                                jj = jj + 1
                                w(i, j, k, l) = w(i, j, k, l) - wdot(jj) * h
                                dwd(i, j, k, l) = aimag(dw(i, j, k, l)) / aimag(h)
                                dwdot(jj) = dwd(i, j, k, l) ! copy values to output
                            end do
                        end do
                    end do
                end do
 
            end do
        end do
 
        alpha = alpha - h * extradot(ialpha)
        beta = beta - h * extradot(ibeta)
        mach = mach - h * extradot(imach)
        machcoef = machcoef - h * extradot(imach)
        machgrid = machgrid - h * extradot(imachgrid)
        pinfdim = pinfdim - h * extradot(ipressure)
        rhoinfdim = rhoinfdim - h * extradot(idensity)
        tinfdim = tinfdim - h * extradot(itemperature)
        pointref(1) = pointref(1) - h * extradot(ipointrefx)
        pointref(2) = pointref(2) - h * extradot(ipointrefy)
        pointref(3) = pointref(3) - h * extradot(ipointrefz)
        rgasdim = rgasdim - h * zero
 
        bcdatavalues = bcdatavalues - bcdatavaluesdot * h
        if (equations == ransequations) then
            call updatexsurf(level)
        end if
 
        fdot = aimag(forces) / aimag(h)
        funcsdot = aimag(funcvalues) / aimag(h)
 
    end subroutine computematrixfreeproductfwdcs
 
#endif
 
end module adjointdebug