MODULE Generic_ViscousAtCellWalls
! ---- The viscous stress terms are calculated on the cell walls
! ---- by means of second order, center difference stencils.
! ---- sigma_ij = mu ( S_ij - 2/3 div.u delta_ij)
! ---- where j = the direction
! ---- and adds this to the flux vector
! ---- note: this requires the same number of ghostcells as the LES
! ---- because the outer-most one is used in calculating sigma
INTERFACE ViscousFluxModification
MODULE PROCEDURE AddViscInterpOneD
MODULE PROCEDURE AddViscInterpTwoD
MODULE PROCEDURE AddViscInterpThreeD
END INTERFACE
CONTAINS
SUBROUTINE AddViscInterpOneD(fxi,ux,vx,direction)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_Transport
USE Generic_AcousticBC
USE Generic_XCellWallAv
USE Generic_XCellWallDx
use cles_interfaces
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(inOUT) :: fxi(ncomps,ixlo:ixhi,1)
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
INTEGER, INTENT(IN) :: direction
DOUBLE PRECISION, ALLOCATABLE :: sigma_1dir(:)
! working arrays
DOUBLE PRECISION, ALLOCATABLE :: cwVisc(:), cwKappa(:), cwDiff(:)
INTEGER :: i, is, iv, iret
DOUBLE PRECISION :: gamma, R, xi(1)
ALLOCATE (sigma_1dir(ixlo:ixhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpOneD/sigma_1dir', iret)
endif
ALLOCATE (cwVisc(ixlo:ixhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpOneD/cwVisc', iret)
endif
ALLOCATE (cwKappa(ixlo:ixhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpOneD/cwKappa', iret)
endif
ALLOCATE (cwDiff(ixlo:ixhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpOneD/cwDiff', iret)
endif
! Viscous term
slct_modev : SELECT CASE (trsp_visc_mode)
CASE (CLES_TRANSP_NONE)
cwVisc = 0.0d0
CASE (CLES_TRANSP_VISC_CST)
cwVisc = mu(1,1,1)
CASE (CLES_TRANSP_VISC_VAR)
CALL XCellWallAv(mu(:,1,1), cwVisc)
END SELECT slct_modev
slct_modek : SELECT CASE (trsp_kappa_mode)
CASE (CLES_TRANSP_NONE)
cwKappa = 0.0d0
CASE (CLES_TRANSP_KAPPA_CST)
cwKappa = kappa(1,1,1)
CASE (CLES_TRANSP_KAPPA_VAR)
CALL XCellWallAv(kappa(:,1,1), cwKappa)
END SELECT slct_modek
! ---- S_11-2/3div.u: 2 du/dx - 2/3 du/dx
! ---- : mu 4/3 du/dx
CALL XCellWallDx(vx, sigma_1dir, nvars, 2)
sigma_1dir = cwVisc*4.0d0/3.0d0*sigma_1dir
CALL AcousticViscousBC(1, sigma_1dir, direction)
! ---- SUBTRACT VISCOUS TERMS FROM THE INTERPOLATED FLUX
! -- fxi(2) - 4/3 du/dx
fxi(2,:,direction) = fxi(2,:,direction) - sigma_1dir(:)
! -- fxi(5) - rho u 4/3 du/dx (use cwDiff as temporary storage)
CALL XCellWallAv(vx, cwDiff, nvars, 2)
fxi(5,:,direction) = fxi(5,:,direction) - cwDiff(:)*sigma_1dir(:)
! ---- Heat Conduction
CALL XCellWallDx(ux, sigma_1dir, ncomps, nvars+1)
! -- kappa* d(T)/dx at the cell wall
sigma_1dir = cwKappa*sigma_1dir
CALL AcousticViscousBC(5, sigma_1dir, direction)
! --- SUBTRACT HEATING TERM FROM INTERPOLATED FLUX
! -- fxi(5) - kappa d(T)/dx
fxi(5,:,direction) = fxi(5,:,direction) - sigma_1dir(:)
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_LEWIS ) then
do i=ixlo, ixhi
! evaluate gamma and R to get cp
call cles_roe(ux(1,i), ux(1,i), ncomps, vx(1,i), vx(1,i), &
nvars, gamma, R, xi, 0)
cwKappa(i) = cwKappa(i)*(gamma-1.0d0)/R/gamma
enddo
endif
! SCALAR DIFFUSION
DO is=1, nscal
iv = is + 5
slct_moded : SELECT CASE (trsp_diff_mode)
CASE (CLES_TRANSP_NONE)
cwDiff = 0.0d0
CASE (CLES_TRANSP_DIFF_CST)
cwDiff = diff(is,1,1,1)
CASE (CLES_TRANSP_DIFF_VAR)
CALL XCellWallAv(diff(:,is,1,1), cwDiff)
CASE (CLES_TRANSP_DIFF_SCHMIDT)
cwDiff = cwVisc * diff(is,1,1,1)
CASE (CLES_TRANSP_DIFF_LEWIS)
cwDiff = cwKappa * diff(is,1,1,1)
END SELECT slct_moded
! --- mu d(vx(is))/dx
CALL XCellWallDx(vx, sigma_1dir, nvars, iv)
sigma_1dir = cwDiff*sigma_1dir
CALL AcousticViscousBC(iv, sigma_1dir, direction)
!-- fxi(is) - rho D d(vx(is))/dx
fxi(iv,:,direction) = fxi(iv,:,direction)-sigma_1dir(:)
END DO
DEALLOCATE (sigma_1dir, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpOneD/sigma_1dir', iret)
endif
DEALLOCATE (cwVisc, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpOneD/cwVisc', iret)
endif
DEALLOCATE (cwKappa, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpOneD/cwKappa', iret)
endif
DEALLOCATE (cwDiff, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpOneD/cwDiff', iret)
endif
END SUBROUTINE AddViscInterpOneD
SUBROUTINE AddViscInterpTwoD(fxi, ux, vx, direction)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_AcousticBC
USE Generic_Transport
USE Generic_XCellWallAv
USE Generic_XCellWallDx
USE Generic_XCellWallDy
USE Generic_YCellWallAv
USE Generic_YCellWallDx
USE Generic_YCellWallDy
use cles_interfaces
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,2)
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
INTEGER, INTENT(IN) :: direction
DOUBLE PRECISION, ALLOCATABLE :: sigma_1dir(:,:)
DOUBLE PRECISION, ALLOCATABLE :: sigma_2dir(:,:)
DOUBLE PRECISION, ALLOCATABLE :: cwVisc(:,:)
DOUBLE PRECISION, ALLOCATABLE :: cwKappa(:,:)
DOUBLE PRECISION, ALLOCATABLE :: cwDiff(:,:)
DOUBLE PRECISION, ALLOCATABLE :: work(:,:)
INTEGER :: is, iv, i, j, iret
DOUBLE PRECISION :: c23, gamma, R, xi(1)
! Viscousity
c23 = 2.0d0/3.0d0
ALLOCATE (sigma_1dir(ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpTwoD/sigma_1dir', iret)
endif
ALLOCATE (sigma_2dir(ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpTwoD/sigma_2dir', iret)
endif
ALLOCATE (cwVisc(ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpTwoD/cwVisc', iret)
endif
ALLOCATE (cwKappa(ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpTwoD/cwKappa', iret)
endif
ALLOCATE (cwDiff(ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpTwoD/cwDiff', iret)
endif
ALLOCATE (work(ixlo:ixhi,iylo:iyhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpTwoD/work', iret)
endif
! Viscous term
slct_modev : SELECT CASE (trsp_visc_mode)
CASE (CLES_TRANSP_NONE)
cwVisc = 0.0d0
CASE (CLES_TRANSP_VISC_CST)
cwVisc = mu(1,1,1)
CASE (CLES_TRANSP_VISC_VAR)
! postpone for direction specific version
END SELECT slct_modev
whichdirection: SELECT CASE(direction)
CASE(1)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
CALL XCellWallAv(mu(:,:,1), cwVisc)
endif
! ---- S_11- 2/3 div.u: 2 du/dx - 2/3 (du/dx + dv/dy)
! ---- : 4/3 du/dx - 2/3 dv/dy
CALL XCellWallDx(vx, sigma_1dir, nvars, 2)
CALL XCellWallDy(vx, work, nvars, 3)
sigma_1dir = cwVisc*(2.0d0*sigma_1dir - work)*c23
! ---- S_21: (dv/dx + du/dy)
CALL XCellWallDy(vx, work, nvars, 2)
CALL XCellWallDx(vx, sigma_2dir, nvars, 3)
sigma_2dir = cwVisc*(sigma_2dir + work)
CALL AcousticViscousBC(1, sigma_1dir, direction)
CALL AcousticViscousBC(2, sigma_2dir, direction)
! --- add to flux
fxi(2,:,:,direction) = fxi(2,:,:,direction) - sigma_1dir(:,:)
fxi(3,:,:,direction) = fxi(3,:,:,direction) - sigma_2dir(:,:)
CALL XCellWallAv(vx, work, nvars, 2)
fxi(5,:,:,direction) = fxi(5,:,:,direction) - work(:,:)*sigma_1dir(:,:)
CALL XCellWallAv(vx, work, nvars, 3)
fxi(5,:,:,direction) = fxi(5,:,:,direction) - work(:,:)*sigma_2dir(:,:)
CASE(2)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
CALL YCellWallAv(mu(:,:,1), cwVisc)
endif
! ---- S_21: (dv/dx + du/dy)
CALL YCellWallDy(vx, work, nvars, 2)
CALL YCellWallDx(vx, sigma_1dir, nvars, 3)
sigma_1dir = cwVisc*(sigma_1dir+work)
! ---- S_22 -2/3 div.u: 2dv/dy-2/3 (du/dx + dv/dy)
! : 4/3 dv/dy-2/3 du/dx
CALL YCellWallDx(vx, work, nvars, 2)
CALL YCellWallDy(vx, sigma_2dir, nvars, 3)
sigma_2dir = cwVisc*(2.0d0*sigma_2dir-work)*c23
CALL AcousticViscousBC(1, sigma_1dir, direction)
CALL AcousticViscousBC(2, sigma_2dir, direction)
! --- add to flux
fxi(2,:,:,direction) = fxi(2,:,:,direction) - sigma_1dir(:,:)
fxi(3,:,:,direction) = fxi(3,:,:,direction) - sigma_2dir(:,:)
CALL YCellWallAv(vx, work, nvars, 2)
fxi(5,:,:,direction) = fxi(5,:,:,direction) - sigma_1dir(:,:)*work(:,:)
CALL YCellWallAv(vx, work, nvars, 3)
fxi(5,:,:,direction) = fxi(5,:,:,direction) - sigma_2dir(:,:)*work(:,:)
END SELECT whichdirection
! Heat Conduction
slct_modek : SELECT CASE (trsp_kappa_mode)
CASE (CLES_TRANSP_NONE)
cwKappa = 0.0d0
CASE (CLES_TRANSP_KAPPA_CST)
cwKappa = kappa(1,1,1)
CASE (CLES_TRANSP_KAPPA_VAR)
! postpone for direction specific version
END SELECT slct_modek
whichdirectionpr: SELECT CASE(direction)
CASE(1)
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
CALL XCellWallAv(kappa(:,:,1), cwKappa)
endif
! kappa d(T)/dx
CALL XCellWallDx(ux, work, ncomps, nvars+1)
CASE(2)
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
CALL YCellWallAv(kappa(:,:,1), cwKappa)
endif
! kappa d(T)/dy
CALL YCellWallDy(ux, work, ncomps, nvars+1)
END SELECT whichdirectionpr
sigma_1dir =cwKappa * work
CALL AcousticViscousBC(5, sigma_1dir, direction)
fxi(5,:,:,direction) = fxi(5,:,:,direction) - sigma_1dir(:,:)
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_LEWIS ) then
do j=iylo, iyhi
do i=ixlo, ixhi
! evaluate gamma and R to get cp
call cles_roe(ux(1,i,j), ux(1,i,j), ncomps, vx(1,i,j), &
vx(1,i,j), nvars, gamma, R, xi, 0)
cwKappa(i,j) = cwKappa(i,j)*(gamma-1.0d0)/R/gamma
enddo
enddo
endif
! Scalar diffusion
DO is=1, nscal
iv = is + 5
slct_moded : SELECT CASE (trsp_diff_mode)
CASE (CLES_TRANSP_NONE)
cwDiff = 0.0d0
CASE (CLES_TRANSP_DIFF_CST)
cwDiff = diff(is,1,1,1)
CASE (CLES_TRANSP_DIFF_VAR)
! postpone for direction specific version
CASE (CLES_TRANSP_DIFF_SCHMIDT)
cwDiff = cwVisc * diff(is,1,1,1)
CASE (CLES_TRANSP_DIFF_LEWIS)
cwDiff = cwKappa * diff(is,1,1,1)
END SELECT slct_moded
whichdirectionsc: SELECT CASE(direction)
CASE(1)
! --- mu dpsi/dx
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
CALL XCellWallAv(diff(:,:,is,1), cwDiff)
endif
CALL XCellWallDx(vx, work, nvars, iv)
CASE(2)
! --- mu dpsi/dy
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
CALL YCellWallAv(diff(:,:,is,1), cwDiff)
endif
CALL YCellWallDy(vx, work, nvars, iv)
END SELECT whichdirectionsc
sigma_1dir = cwDiff * work
CALL AcousticViscousBC(iv, sigma_1dir, direction)
fxi(iv,:,:,direction) = fxi(iv,:,:,direction) - sigma_1dir(:,:)
END DO
DEALLOCATE (sigma_1dir, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpTwoD/sigma_1dir', iret)
endif
DEALLOCATE (sigma_2dir, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpTwoD/sigma_2dir', iret)
endif
DEALLOCATE (cwVisc, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpTwoD/cwVisc', iret)
endif
DEALLOCATE (cwKappa, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpTwoD/cwKappa', iret)
endif
DEALLOCATE (cwDiff, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpTwoD/cwDiff', iret)
endif
DEALLOCATE (work, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpTwoD/work', iret)
endif
END SUBROUTINE AddViscInterpTwoD
SUBROUTINE AddViscInterpThreeD(fxi,ux,vx,direction)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_AcousticBC
USE Generic_Transport
USE Generic_XCellWallAv
USE Generic_XCellWallDx
USE Generic_XCellWallDy
USE Generic_XCellWallDz
USE Generic_YCellWallAv
USE Generic_YCellWallDx
USE Generic_YCellWallDy
USE Generic_YCellWallDz
USE Generic_ZCellWallAv
USE Generic_ZCellWallDx
USE Generic_ZCellWallDy
USE Generic_ZCellWallDz
use cles_interfaces
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT) :: fxi(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi,3)
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER, INTENT(IN) :: direction
DOUBLE PRECISION, ALLOCATABLE :: sigma_1dir(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: sigma_2dir(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: sigma_3dir(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: cwVisc(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: cwKappa(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: cwDiff(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:)
INTEGER :: is, iv, i, j, k, iret
EXTERNAL cles_hook_exist, cles_hook2
INTEGER :: has_hooks, ipar(12), ierr, nn, cles_hook_exist
DOUBLE PRECISION :: c23, gamma, R, xi(1)
call cleslog_log_enter('AddViscInterpThreeD')
has_hooks = cles_hook_exist(CLES_HOOK_VISCOUS)
c23 = 2.0d0/3.0d0
ALLOCATE (sigma_1dir(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/sigma_1dir', iret)
endif
ALLOCATE (sigma_2dir(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/sigma_2dir', iret)
endif
ALLOCATE (sigma_3dir(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/sigma_3dir', iret)
endif
ALLOCATE (cwVisc(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/cwVisc', iret)
endif
ALLOCATE (cwKappa(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/cwKappa', iret)
endif
ALLOCATE (cwDiff(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/cwDiff', iret)
endif
ALLOCATE (work(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AddViscInterpThreeD/work', iret)
endif
! Viscous term
slct_modev : SELECT CASE (trsp_visc_mode)
CASE (CLES_TRANSP_NONE)
cwVisc = 0.0d0
CASE (CLES_TRANSP_VISC_CST)
cwVisc = mu(1,1,1)
CASE (CLES_TRANSP_VISC_VAR)
! postpone for direction specific version
END SELECT slct_modev
ipar(1) = nx
ipar(2) = ny
ipar(3) = nz
ipar(4) = ixlo
ipar(5) = ixhi
ipar(6) = iylo
ipar(7) = iyhi
ipar(8) = izlo
ipar(9) = izhi
nn = mx*my*mz
ipar(10) = direction
whichdirection: SELECT CASE(direction)
CASE (1)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
CALL XCellWallAv(mu, cwVisc)
endif
! ---- S_11-2/3div.u: 2du/dx -2/3 (du/dx + dv/dy + dw/dz)
! ---- S_11-2/3div.u: 4/3 du/dx -2/3 (dv/dy + dw/dz)
CALL XCellWallDy(vx, sigma_1dir, nvars, 3)
CALL XCellWallDx(vx, sigma_2dir, nvars, 3)
CALL XCellWallDz(vx, work, nvars, 4)
CALL XCellWallDx(vx, sigma_3dir, nvars, 4)
sigma_1dir = sigma_1dir + work
CALL XCellWallDx(vx, work, nvars, 2)
sigma_1dir = cwVisc * (2.0d0*work - sigma_1dir)*c23
! ---- S_21: (dv/dx + du/dy)
CALL XCellWallDy(vx, work, nvars, 2)
sigma_2dir = cwVisc * (sigma_2dir + work)
! ---- S_31: (dw/dx + du/dz)
CALL XCellWallDz(vx, work, nvars, 2)
sigma_3dir = cwVisc * (sigma_3dir + work)
CALL AcousticViscousBC(1, sigma_1dir, direction)
CALL AcousticViscousBC(2, sigma_2dir, direction)
CALL AcousticViscousBC(3, sigma_3dir, direction)
if ( has_hooks .eq. CLES_TRUE ) then
ipar(11) = 1
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_1dir, nn, ierr)
ipar(11) = 2
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_2dir, nn, ierr)
ipar(11) = 3
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_3dir, nn, ierr)
endif
! --- add to flux
fxi(2,:,:,:,direction) = fxi(2,:,:,:,direction) - sigma_1dir(:,:,:)
fxi(3,:,:,:,direction) = fxi(3,:,:,:,direction) - sigma_2dir(:,:,:)
fxi(4,:,:,:,direction) = fxi(4,:,:,:,direction) - sigma_3dir(:,:,:)
CALL XCellWallAv(vx, work, nvars, 2)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_1dir * work
CALL XCellWallAv(vx, work, nvars, 3)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_2dir * work
CALL XCellWallAv(vx, work, nvars, 4)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_3dir * work
CASE (2)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
CALL YCellWallAv(mu, cwVisc)
endif
! ---- S_12: (du/dy + dv/dx)
CALL YCellWallDy(vx, sigma_1dir, nvars, 2)
CALL YCellWallDx(vx, sigma_2dir, nvars, 2)
CALL YCellWallDz(vx, work, nvars, 4)
sigma_2dir = sigma_2dir + work
CALL YCellWallDy(vx, sigma_3dir, nvars, 4)
CALL YCellWallDx(vx, work, nvars, 3)
sigma_1dir = cwVisc*(sigma_1dir + work)
! ---- S_22-2/3div.u: 2(dv/dy) -2/3 (du/dx + dv/dy + dw/dz)
! ---- : 4/3 (dv/dy) -2/3 (du/dx + dw/dz)
CALL YCellWallDy(vx, work, nvars, 3)
sigma_2dir = cwVisc*(2.0d0*work - sigma_2dir)*c23
! ---- S_32: (dw/dy + dv/dz)
CALL YCellWallDz(vx, work, nvars, 3)
sigma_3dir = cwVisc*(sigma_3dir + work)
CALL AcousticViscousBC(1, sigma_1dir, direction)
CALL AcousticViscousBC(2, sigma_2dir, direction)
CALL AcousticViscousBC(3, sigma_3dir, direction)
if ( has_hooks .eq. CLES_TRUE ) then
ipar(11) = 1
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_1dir, nn, ierr)
ipar(11) = 2
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_2dir, nn, ierr)
ipar(11) = 3
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_3dir, nn, ierr)
endif
! --- add to flux
fxi(2,:,:,:,direction) = fxi(2,:,:,:,direction) - sigma_1dir(:,:,:)
fxi(3,:,:,:,direction) = fxi(3,:,:,:,direction) - sigma_2dir(:,:,:)
fxi(4,:,:,:,direction) = fxi(4,:,:,:,direction) - sigma_3dir(:,:,:)
CALL YCellWallAv(vx, work, nvars, 2)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_1dir * work
CALL YCellWallAv(vx, work, nvars, 3)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_2dir * work
CALL YCellWallAv(vx, work, nvars, 4)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_3dir * work
CASE (3)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
CALL ZCellWallAv(mu, cwVisc)
endif
! ---- S_13: (du/dz + dw/dx)
CALL ZCellWallDz(vx, sigma_1dir, nvars, 2)
CALL ZCellWallDx(vx, sigma_3dir, nvars, 2)
CALL ZCellWallDy(vx, work, nvars, 3)
sigma_3dir = sigma_3dir + work
CALL ZCellWallDz(vx, sigma_2dir, nvars, 3)
CALL ZCellWallDx(vx, work, nvars, 4)
sigma_1dir = cwVisc*(sigma_1dir + work)
! ---- S_23: (dw/dy + dv/dz)
CALL ZCellWallDy(vx, work, nvars, 4)
sigma_2dir = cwVisc*(sigma_2dir + work)
! ---- S_33-2/3div.u: 2(dw/dz) -2/3 (du/dx + dv/dy + dw/dz)
! ---- : 4/3(dw/dz) -2/3 (du/dx + dv/dy)
CALL ZCellWallDz(vx, work, nvars, 4)
sigma_3dir = cwVisc*(2.0d0*work - sigma_3dir)*c23
CALL AcousticViscousBC(1, sigma_1dir, direction)
CALL AcousticViscousBC(2, sigma_2dir, direction)
CALL AcousticViscousBC(3, sigma_3dir, direction)
if ( has_hooks .eq. CLES_TRUE ) then
ipar(11) = 1
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_1dir, nn, ierr)
ipar(11) = 2
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_2dir, nn, ierr)
ipar(11) = 3
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_3dir, nn, ierr)
endif
! --- add to flux
fxi(2,:,:,:,direction) = fxi(2,:,:,:,direction) - sigma_1dir(:,:,:)
fxi(3,:,:,:,direction) = fxi(3,:,:,:,direction) - sigma_2dir(:,:,:)
fxi(4,:,:,:,direction) = fxi(4,:,:,:,direction) - sigma_3dir(:,:,:)
CALL ZCellWallAv(vx, work, nvars, 2)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_1dir * work
CALL ZCellWallAv(vx, work, nvars, 3)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_2dir * work
CALL ZCellWallAv(vx, work, nvars, 4)
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_3dir * work
END SELECT whichdirection
! Heat Conduction
slct_modek : SELECT CASE (trsp_kappa_mode)
CASE (CLES_TRANSP_NONE)
cwKappa = 0.0d0
CASE (CLES_TRANSP_KAPPA_CST)
cwKappa = kappa(1,1,1)
CASE (CLES_TRANSP_KAPPA_VAR)
! postpone for direction specific version
END SELECT slct_modek
whichdirectionpr: SELECT CASE(direction)
CASE(1)
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
CALL XCellWallAv(kappa, cwKappa)
endif
! kappa d(T)/dx
CALL XCellWallDx(ux, work, ncomps, nvars+1)
CASE(2)
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
CALL YCellWallAv(kappa, cwKappa)
endif
! kappa d(T)/dy
CALL YCellWallDy(ux, work, ncomps, nvars+1)
CASE(3)
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
CALL ZCellWallAv(kappa, cwKappa)
endif
! kappa d(T)/dy
CALL ZCellWallDz(ux, work, ncomps, nvars+1)
END SELECT whichdirectionpr
sigma_1dir = cwKappa * work
CALL AcousticViscousBC(5, sigma_1dir, direction)
if ( has_hooks .eq. CLES_TRUE ) then
ipar(11) = 4
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_1dir, nn, ierr)
endif
fxi(5,:,:,:,direction) = fxi(5,:,:,:,direction) - sigma_1dir(:,:,:)
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_LEWIS ) then
do k=izlo, izhi
do j=iylo, iyhi
do i=ixlo, ixhi
! evaluate gamma and R to get cp
call cles_roe(ux(1,i,j,k), ux(1,i,j,k), ncomps, vx(1,i,j,k), &
vx(1,i,j,k), nvars, gamma, R, xi, 0)
cwKappa(i,j,k) = cwKappa(i,j,k)*(gamma-1.0d0)/R/gamma
enddo
enddo
enddo
endif
! Scalar diffusion
DO is=1, nscal
iv = is + 5
slct_moded : SELECT CASE (trsp_diff_mode)
CASE (CLES_TRANSP_NONE)
cwDiff = 0.0d0
CASE (CLES_TRANSP_DIFF_CST)
cwDiff = diff(is,1,1,1)
CASE (CLES_TRANSP_DIFF_VAR)
! postpone
CASE (CLES_TRANSP_DIFF_SCHMIDT)
cwDiff = cwVisc * diff(is,1,1,1)
CASE (CLES_TRANSP_DIFF_LEWIS)
cwDiff = cwKappa * diff(is,1,1,1)
END SELECT slct_moded
whichdirectionsc: SELECT CASE(direction)
CASE(1)
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
CALL XCellWallAv(diff(:,:,:,is), cwDiff)
endif
! mu d(psi)/dx
CALL XCellWallDx(vx, work, nvars, iv)
CASE(2)
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
CALL YCellWallAv(diff(:,:,:,is), cwDiff)
endif
! mu d(psi)/dy
CALL YCellWallDy(vx, work, nvars, iv)
CASE(3)
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
CALL ZCellWallAv(diff(:,:,:,is), cwDiff)
endif
! mu d(psi)/dz
CALL ZCellWallDz(vx, work, nvars, iv)
END SELECT whichdirectionsc
sigma_1dir = cwDiff * work
CALL AcousticViscousBC(iv, sigma_1dir, direction)
if ( has_hooks .eq. CLES_TRUE ) then
ipar(11) = iv
call cles_hook2(CLES_HOOK_VISCOUS, ipar, 11, sigma_1dir, nn, ierr)
endif
fxi(iv,:,:,:,direction) = fxi(iv,:,:,:,direction) - sigma_1dir(:,:,:)
END DO
DEALLOCATE (sigma_1dir, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/sigma_1dir', iret)
endif
DEALLOCATE (sigma_2dir, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/sigma_2dir', iret)
endif
DEALLOCATE (sigma_3dir, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/sigma_3dir', iret)
endif
DEALLOCATE (cwVisc, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/cwVisc', iret)
endif
DEALLOCATE (cwKappa, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/cwKappa', iret)
endif
DEALLOCATE (cwDiff, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/cwDiff', iret)
endif
DEALLOCATE (work, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'AddViscInterpThreeD/work', iret)
endif
call cleslog_log_exit('AddViscInterpThreeD')
END SUBROUTINE AddViscInterpThreeD
END Module Generic_ViscousAtCellWalls