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