smoothers.F90

Go to the documentation of this file.

module smoothers
 
contains
    subroutine rungekuttasmoother
        !
        !       RungeKuttaSmoother performs one multi-stage runge kutta
        !       explicit time step for the current multigrid level. On
        !       entrance it is assumed that the residual and time step are
        !       already computed. On exit the solution in the halo's contain
        !       the latest values. However, the residual corresponding to
        !       these values is not computed.
        !
        use constants
        use blockpointers, only: w, p, wn, pn, il, jl, kl, ndom
        use flowvarrefstate, only: nwf
        use inputiteration, only: nrkstages
        use inputtimespectral, only: ntimeintervalsspectral
        use iteration, only: currentlevel, rkstage
        use utils, only: setpointers
        use residuals, only: initres, residual, residualaveraging, sourceterms
        use bcroutines, only: applyallbc
        implicit none
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, nn, i, j, k, l
 
        ! Store the variables of the zeroth runge kutta stage.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
            domains: do nn = 1, ndom
 
                ! Set the pointers to this block.
 
                call setpointers(nn, currentlevel, sps)
 
                ! The variables stored in w.
 
                do l = 1, nwf
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
                                wn(i, j, k, l) = w(i, j, k, l)
                            end do
                        end do
                    end do
                end do
 
                ! And the pressure.
 
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
                            pn(i, j, k) = p(i, j, k)
                        end do
                    end do
                end do
 
            end do domains
        end do spectralloop
 
        ! Loop over all but the last Runge Kutta stages. Note that the
        ! counter variable rkStage is defined in the module iteration.
 
        do rkstage = 1, (nrkstages - 1)
 
            ! Execute a Runge Kutta stage and exchange the externals.
 
            call executerkstage
 
            ! Compute the residuals for the next stage.
 
            call initres(1_inttype, nwf)
            call sourceterms()
            call residual
 
        end do
 
        ! Execute the last RK stage. Set rkStage to nRKStages, for
        ! clarity; after the previous loop rkStage == nRKStages.
 
        rkstage = nrkstages
 
        call executerkstage
 
    end subroutine rungekuttasmoother
 
    !      ==================================================================
 
    subroutine executerkstage
        !
        !       executeRkStage executes one runge kutta stage. The stage
        !       number, rkStage, is defined in the local module iteration.
        !
        use blockpointers
        use constants
        use flowvarrefstate
        use inputiteration
        use inputphysics
        use inputtimespectral
        use inputunsteady
        use iteration
        use inputdiscretization
        use haloexchange, only: whalo1, whalo2
        use utils, only: setpointers
        use flowutils, only: computeetotblock, computelamviscosity
        use turbutils, only: computeeddyviscosity
        use residuals, only: residualaveraging
        use bcroutines, only: applyallbc
 
        implicit none
        !
        !      Local parameter.
        !
        real(kind=realtype), parameter :: fivethird = five * third
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, nn, i, j, k, l
 
        real(kind=realtype) :: tmp, unsteadyimpl, mult
        real(kind=realtype) :: dt, currentcfl, gm1, gm53
        real(kind=realtype) :: v2, ovr, dp, factk, ru, rv, rw
 
        logical :: secondHalo, smoothResidual, correctForK
 
        ! Set the value of secondHalo and the current cfl number,
        ! depending on the situation. On the finest grid in the mg cycle
        ! the second halo is computed, otherwise not.
 
        if (currentlevel <= groundlevel) then
            secondhalo = .true.
        else
            secondhalo = .false.
        end if
 
        currentcfl = cflcoarse
        if (currentlevel == 1) then
            currentcfl = cfl
        end if
 
        ! Determine whether or not residual averaging must be applied.
 
        if (resaveraging == noresaveraging) then
            smoothresidual = .false.
        else if (resaveraging == alwaysresaveraging) then
            smoothresidual = .true.
        else if (mod(rkstage, 2_inttype) == 1) then
            smoothresidual = .true.
        else
            smoothresidual = .false.
        end if
 
        ! Determine whether or not the total energy must be corrected
        ! for the presence of the turbulent kinetic energy.
 
        if (kpresent) then
            if ((currentlevel <= groundlevel)) then
                correctfork = .true.
            else
                correctfork = .false.
            end if
        else
            correctfork = .false.
        end if
        !
        !       Compute the updates of the conservative variables.
        !
        ! Loop over the local number of blocks.
 
        domainsupdate: do nn = 1, ndom
 
            ! Determine the equation mode solved.
 
            select case (equationmode)
 
            case (steady, timespectral)
 
                ! Steady equations, including time spectral. Everything is
                ! solved explicitly. Store the cfl number times the RK
                ! coefficient in tmp.
 
                tmp = currentcfl * etark(rkstage)
 
                ! Loop over the number of spectral solutions. Note that
                ! for the steady mode this value is 1.
 
                spectralsteady: do sps = 1, ntimeintervalsspectral
 
                    ! Set the pointers to this block.
 
                    call setpointers(nn, currentlevel, sps)
 
                    ! Loop over the owned cells of this block.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
 
                                ! Determine the local time step (multiplied by the
                                ! current rk coefficient).
                                if (lowspeedpreconditioner) then
                                    dt = 0.8 * tmp * dtl(i, j, k)
                                else
                                    dt = tmp * dtl(i, j, k)
                                end if
 
                                ! Compute the updates of the flow field variables.
 
                                dw(i, j, k, irho) = dw(i, j, k, irho) * dt
                                dw(i, j, k, imx) = dw(i, j, k, imx) * dt
                                dw(i, j, k, imy) = dw(i, j, k, imy) * dt
                                dw(i, j, k, imz) = dw(i, j, k, imz) * dt
                                dw(i, j, k, irhoe) = dw(i, j, k, irhoe) * dt
 
                            end do
                        end do
                    end do
 
                end do spectralsteady
 
                !=============================================================
 
            case (unsteady)
 
                ! Unsteady equations are solved via the dual time
                ! stepping technique. The leading term of the time
                ! integrator must be treated implicitly for stability
                ! reasons. This leads to a different multiplication
                ! factor of the residual compared to the steady case.
                ! Compute this additional term and store the cfl number
                ! times the rk coefficient in tmp.
 
                unsteadyimpl = coeftime(0) * timeref / deltat
                tmp = currentcfl * etark(rkstage)
 
                ! Loop over the number of spectral modes, although this is
                ! always 1 for the unsteady mode. The loop is executed for
                ! consistency reasons.
 
                spectralunsteady: do sps = 1, ntimeintervalsspectral
 
                    ! Set the pointers to this block.
 
                    call setpointers(nn, currentlevel, sps)
 
                    ! Determine the updates of the flow field variables.
                    ! Owned cells only. The result is stored in dw.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
 
                                ! Determine the local time step (multiplied by the
                                ! current rk coefficient) and the multiplication
                                ! factor for the residuals.
 
                                dt = tmp * dtl(i, j, k)
                                mult = dt / (dt * unsteadyimpl * vol(i, j, k) + one)
 
                                ! Compute the updates of the flow field variables.
 
                                dw(i, j, k, irho) = dw(i, j, k, irho) * mult
                                dw(i, j, k, imx) = dw(i, j, k, imx) * mult
                                dw(i, j, k, imy) = dw(i, j, k, imy) * mult
                                dw(i, j, k, imz) = dw(i, j, k, imz) * mult
                                dw(i, j, k, irhoe) = dw(i, j, k, irhoe) * mult
 
                            end do
                        end do
                    end do
 
                end do spectralunsteady
 
            end select
 
        end do domainsupdate
        !
        !       Compute the new state vector.
        !
        ! Loop over the number of spectral solutions and local blocks.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
            domainsstate: do nn = 1, ndom
 
                ! Set the pointers to this block.
 
                call setpointers(nn, currentlevel, sps)
 
                ! Possibility to smooth the updates.
 
                if (smoothresidual) then
                    call residualaveraging
                end if
                ! Flow variables.
 
                factk = zero
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
 
                            ! Store gamma -1 and gamma - 5/3 a bit easier.
 
                            gm1 = gamma(i, j, k) - one
                            gm53 = gamma(i, j, k) - fivethird
 
                            ! Compute the pressure update from the conservative
                            ! updates. The expression used below is valid even if
                            ! cp is not constant. For the calorically perfect case,
                            ! cp is constant, it can be simplified to the usual
                            ! expression, but not for variable cp.
 
                            ovr = one / w(i, j, k, irho)
                            v2 = w(i, j, k, ivx)**2 + w(i, j, k, ivy)**2 + w(i, j, k, ivz)**2
                            if (correctfork) factk = gm53 * w(i, j, k, itu1)
 
                            dp = (ovr * p(i, j, k) + factk &
                                  - gm1 * (ovr * w(i, j, k, irhoe) - v2)) * dw(i, j, k, irho) &
                                 + gm1 * (dw(i, j, k, irhoe) - w(i, j, k, ivx) * dw(i, j, k, imx) &
                                          - w(i, j, k, ivy) * dw(i, j, k, imy) &
                                          - w(i, j, k, ivz) * dw(i, j, k, imz))
 
                            ! Compute the density. Correct for negative values.
 
                            w(i, j, k, irho) = wn(i, j, k, irho) - dw(i, j, k, irho)
                            w(i, j, k, irho) = max(w(i, j, k, irho), 1.e-4_realtype * rhoinf)
 
                            ! Compute the velocities.
 
                            ru = wn(i, j, k, irho) * wn(i, j, k, ivx) - dw(i, j, k, imx)
                            rv = wn(i, j, k, irho) * wn(i, j, k, ivy) - dw(i, j, k, imy)
                            rw = wn(i, j, k, irho) * wn(i, j, k, ivz) - dw(i, j, k, imz)
 
                            ovr = one / w(i, j, k, irho)
                            w(i, j, k, ivx) = ovr * ru
                            w(i, j, k, ivy) = ovr * rv
                            w(i, j, k, ivz) = ovr * rw
 
                            ! Compute the pressure. Correct for negative values.
 
                            p(i, j, k) = pn(i, j, k) - dp
                            p(i, j, k) = max(p(i, j, k), 1.e-4_realtype * pinfcorr)
 
                        end do
                    end do
                end do
 
                ! Compute the total energy and possibly the laminar and eddy
                ! viscosity in the owned cells.
 
                call computeetotblock(2_inttype, il, 2_inttype, jl, &
                                      2_inttype, kl, correctfork)
                call computelamviscosity(.false.)
                call computeeddyviscosity(.false.)
 
            end do domainsstate
        end do spectralloop
 
        ! Exchange the pressure if the pressure must be exchanged early.
        ! Only the first halo's are needed, thus whalo1 is called.
        ! Only on the fine grid.
 
        if (exchangepressureearly .and. currentlevel <= groundlevel) &
            call whalo1(currentlevel, 1_inttype, 0_inttype, .true., &
                        .false., .false.)
 
        ! Apply all boundary conditions to all blocks on this level.
 
        call applyallbc(secondhalo)
 
        ! Exchange the solution. Either whalo1 or whalo2
        ! must be called.
 
        if (secondhalo) then
            call whalo2(currentlevel, 1_inttype, nwf, .true., &
                        .true., .true.)
        else
            call whalo1(currentlevel, 1_inttype, nwf, .true., &
                        .true., .true.)
        end if
 
    end subroutine executerkstage
    subroutine dadismoother
        !
        !       RungeKuttaSmoother performs one multi-stage runge kutta
        !       explicit time step for the current multigrid level. On
        !       entrance it is assumed that the residual and time step are
        !       already computed. On exit the solution in the halo's contain
        !       the latest values. However, the residual corresponding to
        !       these values is not computed.
        !
        use blockpointers
        use flowvarrefstate
        use inputiteration
        use inputtimespectral
        use iteration
        use residuals, only: initres, residual, sourceterms
        implicit none
 
        if (groundlevel == 1) then
            do subit = 1, nsubiterations - 1
 
                ! Execute a DADI step and exchange the externals.
 
                call executedadistep
 
                ! Compute the residuals for the next stage.
                call initres(1_inttype, nwf)
                call sourceterms()
                call residual
 
            end do
 
            ! Set Subit to nSubiterations, for clarity;
            subit = nsubiterations
        end if
 
        ! Execute the last subiteration.
        call executedadistep
 
    end subroutine dadismoother
 
    !      ==================================================================
 
    subroutine executedadistep
        !
        !       executeDADIStep executes one DADI step.
        !
        use blockpointers
        use constants
        use flowvarrefstate
        use inputiteration
        use inputphysics
        use inputtimespectral
        use inputunsteady
        use iteration
        use utils, only: getcorrectfork, setpointers
        use haloexchange, only: whalo1, whalo2
        use flowutils, only: computeetotblock, computelamviscosity
        use turbutils, only: computeeddyviscosity
        use residuals, only: residualaveraging, computedwdadi
        use bcroutines, only: applyallbc
        implicit none
        !
        !      Local parameter.
        !
        real(kind=realtype), parameter :: fivethird = five * third
        !
        !      Local variables.
        !
        integer(kind=intType) :: sps, nn, i, j, k, l
 
        real(kind=realtype) :: unsteadyimpl, mult
        real(kind=realtype) :: dt, currentcfl, gm1, gm53
        real(kind=realtype) :: v2, ovr, dp, factk, ru, rv, rw
 
        logical :: secondHalo, smoothResidual, correctForK
 
        ! Set the value of secondHalo and the current cfl number,
        ! depending on the situation. On the finest grid in the mg cycle
        ! the second halo is computed, otherwise not.
 
        if (currentlevel <= groundlevel) then
            secondhalo = .true.
        else
            secondhalo = .false.
        end if
 
        currentcfl = cflcoarse
        if (currentlevel == 1) then
            currentcfl = cfl
        end if
 
        ! Determine whether or not residual averaging must be applied.
 
        if (resaveraging == noresaveraging) then
            smoothresidual = .false.
        else if (resaveraging == alwaysresaveraging) then
            smoothresidual = .true.
        else if (mod(rkstage, 2_inttype) == 1) then
            smoothresidual = .true.
        else
            smoothresidual = .false.
        end if
 
        ! Determine whether or not the total energy must be corrected
        ! for the presence of the turbulent kinetic energy.
        correctfork = getcorrectfork()
 
        !
        !       Compute the updates of the conservative variables.
        !
        ! Loop over the local number of blocks.
 
        domainsupdate: do nn = 1, ndom
 
            ! Determine the equation mode solved.
 
            select case (equationmode)
 
            case (steady, timespectral)
 
                ! Loop over the number of spectral solutions. Note that
                ! for the steady mode this value is 1.
 
                spectralsteady: do sps = 1, ntimeintervalsspectral
 
                    ! Set the pointers to this block.
 
                    call setpointers(nn, currentlevel, sps)
 
                    ! Loop over the owned cells of this block.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
 
                                ! Determine the local time step
 
                                dt = -currentcfl * dtl(i, j, k) * vol(i, j, k)
 
                                ! Compute the updates of the flow field variables.
 
                                dw(i, j, k, irho) = dw(i, j, k, irho) * dt
                                dw(i, j, k, imx) = dw(i, j, k, imx) * dt
                                dw(i, j, k, imy) = dw(i, j, k, imy) * dt
                                dw(i, j, k, imz) = dw(i, j, k, imz) * dt
                                dw(i, j, k, irhoe) = dw(i, j, k, irhoe) * dt
 
                            end do
                        end do
                    end do
 
                    call computedwdadi
 
                end do spectralsteady
 
                !=============================================================
 
            case (unsteady)
 
                ! Unsteady equations are solved via the dual time
                ! stepping technique. The leading term of the time
                ! integrator must be treated implicitly for stability
                ! reasons. This leads to a different multiplication
                ! factor of the residual compared to the steady case.
 
                unsteadyimpl = coeftime(0) * timeref / deltat
 
                ! Loop over the number of spectral modes, although this is
                ! always 1 for the unsteady mode. The loop is executed for
                ! consistency reasons.
 
                spectralunsteady: do sps = 1, ntimeintervalsspectral
 
                    ! Set the pointers to this block.
 
                    call setpointers(nn, currentlevel, sps)
 
                    ! Determine the updates of the flow field variables.
                    ! Owned cells only. The result is stored in dw.
 
                    do k = 2, kl
                        do j = 2, jl
                            do i = 2, il
 
                                ! Determine the local time step
 
                                dt = currentcfl * dtl(i, j, k)
                                mult = dt / (dt * unsteadyimpl * vol(i, j, k) + one)
                                mult = -mult * vol(i, j, k)
 
                                dw(i, j, k, irho) = dw(i, j, k, irho) * mult
                                dw(i, j, k, imx) = dw(i, j, k, imx) * mult
                                dw(i, j, k, imy) = dw(i, j, k, imy) * mult
                                dw(i, j, k, imz) = dw(i, j, k, imz) * mult
                                dw(i, j, k, irhoe) = dw(i, j, k, irhoe) * mult
 
                            end do
                        end do
                    end do
 
                    call computedwdadi
 
                end do spectralunsteady
 
            end select
 
        end do domainsupdate
        !
        !       Compute the new state vector.
        !
        ! Loop over the number of spectral solutions and local blocks.
 
        spectralloop: do sps = 1, ntimeintervalsspectral
            domainsstate: do nn = 1, ndom
 
                ! Set the pointers to this block.
 
                call setpointers(nn, currentlevel, sps)
 
                ! Possibility to smooth the updates.
 
                if (smoothresidual) call residualaveraging
 
                ! Flow variables.
 
                factk = zero
                do k = 2, kl
                    do j = 2, jl
                        do i = 2, il
 
                            ! Store gamma -1 and gamma - 5/3 a bit easier.
 
                            gm1 = gamma(i, j, k) - one
                            gm53 = gamma(i, j, k) - fivethird
 
                            ! Compute the pressure update from the conservative
                            ! updates. The expression used below is valid even if
                            ! cp is not constant. For the calorically perfect case,
                            ! cp is constant, it can be simplified to the usual
                            ! expression, but not for variable cp.
 
                            ovr = one / w(i, j, k, irho)
                            v2 = w(i, j, k, ivx)**2 + w(i, j, k, ivy)**2 + w(i, j, k, ivz)**2
                            if (correctfork) factk = gm53 * w(i, j, k, itu1)
 
                            dp = (ovr * p(i, j, k) + factk &
                                  - gm1 * (ovr * w(i, j, k, irhoe) - v2)) * dw(i, j, k, irho) &
                                 + gm1 * (dw(i, j, k, irhoe) - w(i, j, k, ivx) * dw(i, j, k, imx) &
                                          - w(i, j, k, ivy) * dw(i, j, k, imy) &
                                          - w(i, j, k, ivz) * dw(i, j, k, imz))
 
                            ! Compute the velocities.
 
                            ru = w(i, j, k, irho) * w(i, j, k, ivx) - dw(i, j, k, imx)
                            rv = w(i, j, k, irho) * w(i, j, k, ivy) - dw(i, j, k, imy)
                            rw = w(i, j, k, irho) * w(i, j, k, ivz) - dw(i, j, k, imz)
 
                            ! Compute the density. Correct for negative values.
 
                            w(i, j, k, irho) = w(i, j, k, irho) - dw(i, j, k, irho)
                            w(i, j, k, irho) = max(w(i, j, k, irho), 1.e-4_realtype * rhoinf)
 
                            ovr = one / w(i, j, k, irho)
                            w(i, j, k, ivx) = ovr * ru
                            w(i, j, k, ivy) = ovr * rv
                            w(i, j, k, ivz) = ovr * rw
 
                            ! Compute the pressure. Correct for negative values.
 
                            p(i, j, k) = p(i, j, k) - dp
                            p(i, j, k) = max(p(i, j, k), 1.e-4_realtype * pinfcorr)
 
                        end do
                    end do
                end do
 
                ! Compute the total energy and possibly the laminar and eddy
                ! viscosity in the owned cells.
 
                call computeetotblock(2_inttype, il, 2_inttype, jl, &
                                      2_inttype, kl, correctfork)
                call computelamviscosity(.false.)
                call computeeddyviscosity(.false.)
 
            end do domainsstate
        end do spectralloop
 
        ! Exchange the pressure if the pressure must be exchanged early.
        ! Only the first halo's are needed, thus whalo1 is called.
        ! Only on the fine grid.
 
        if (exchangepressureearly .and. currentlevel <= groundlevel) &
            call whalo1(currentlevel, 1_inttype, 0_inttype, .true., &
                        .false., .false.)
 
        ! Apply all boundary conditions to all blocks on this level.
 
        call applyallbc(secondhalo)
 
        ! Exchange the solution. Either whalo1 or whalo2
        ! must be called.
 
        if (secondhalo) then
            call whalo2(currentlevel, 1_inttype, nwf, .true., &
                        .true., .true.)
        else
            call whalo1(currentlevel, 1_inttype, nwf, .true., &
                        .true., .true.)
        end if
 
    end subroutine executedadistep
 
end module smoothers