! $Id$
!
! This modules implements viscous heating and diffusion terms
! here smagorinsky viscosity
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: lviscosity = .true.
!
! MVAR CONTRIBUTION 0
! MAUX CONTRIBUTION 1
!
!***************************************************************
module Viscosity
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Density
use Messages
!
implicit none
!
character (len=labellen) :: ivisc='smagorinsky'
integer :: idiag_dtnu=0
real :: maxeffectivenu,nu_mol,C_smag=0.0
logical :: lvisc_first=.false.
!
integer :: idiag_nu_LES=0 ! DIAG_DOC: Mean value of Smagorinsky viscosity
!
namelist /viscosity_run_pars/ nu, lvisc_first,ivisc,c_smag
!
contains
!***********************************************************************
subroutine register_viscosity()
!
! 19-nov-02/tony: coded
! 16-july-04/nils: adapted from visc_shock
!
use Cdata
use FArrayManager
use Sub
!
lvisc_smagorinsky = .true.
!
call farray_register_auxiliary('smagorinsky',ismagorinsky,communicated=.true.)
!
! identify version number
!
if (lroot) call svn_id( &
"$Id$")
!
! Writing files for use with IDL
!
if (naux < maux) aux_var(aux_count)=',smagorinsky $'
if (naux == maux) aux_var(aux_count)=',smagorinsky'
aux_count=aux_count+1
if (lroot) write(15,*) 'smagorinsky = fltarr(mx,my,mz)*one'
!
endsubroutine register_viscosity
!***********************************************************************
subroutine initialize_viscosity()
!
! 20-nov-02/tony: coded
! 16-july-04/nils: adapted from visc_shock
!
use Cdata
!
if (headtt.and.lroot) print*,'viscosity: nu=',nu
if (headtt.and.lroot) print*,'viscosity: c_smag=',c_smag
!
lneed_sij=.true.
if (.not.ldensity) &
call warning('initialize_viscosity',"ldensity better be .true. for ivisc='smagorinsky*'")
!
endsubroutine initialize_viscosity
!*******************************************************************
subroutine rprint_viscosity(lreset,lwrite)
!
! Writes ismagorinsky to index.pro file
!
! 16-july-04/nils: adapted from visc_shock
!
use Cdata
use Diagnostics
use Sub
!
logical :: lreset
logical, optional :: lwrite
integer :: iname
!
! reset everything in case of reset
! (this needs to be consistent with what is defined above!)
!
if (lreset) then
idiag_nu_LES=0
endif
!
! iname runs through all possible names that may be listed in print.in
!
if (lroot.and.ip<14) print*,'rprint_ionization: run through parse list'
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'nu_LES',idiag_nu_LES)
enddo
!
! write column where which viscosity variable is stored
!
if (present(lwrite)) then
if (lwrite) then
!
endif
endif
!
call keep_compiler_quiet(lreset)
!
endsubroutine rprint_viscosity
!***********************************************************************
subroutine pencil_criteria_viscosity()
!
! All pencils that the Viscosity module depends on are specified here.
!
! 21-11-04/anders: coded
!
endsubroutine pencil_criteria_viscosity
!***********************************************************************
subroutine pencil_interdep_viscosity(lpencil_in)
!
! Interdependency among pencils from the Viscosity module is specified here.
!
! 21-11-04/anders: coded
!
use Cdata
!
logical, dimension (npencils) :: lpencil_in
!
call keep_compiler_quiet(lpencil_in)
!
endsubroutine pencil_interdep_viscosity
!***********************************************************************
subroutine calc_pencils_viscosity(f,p)
!
! Calculate Viscosity pencils.
! Most basic pencils should come first, as others may depend on them.
!
! 21-11-04/anders: coded
!
use Cdata
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
!
intent(in) :: f,p
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(p)
!
endsubroutine calc_pencils_viscosity
!***********************************************************************
subroutine calc_viscosity(f)
!
! calculate nu_smag for full domain
!
! 16-july-04/nils: coded
!
use Diagnostics
use Mpicomm
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,2) :: tmp
real, dimension (nx,ny,nz) :: sij2
real, dimension (nx) :: nu_smag
integer :: i,j,ncount,mcount
!
!
!
sij2=0
do i=1,2
do j=i+1,3
!
! find uij and uji (stored in tmp in order to save space)
!
call der_2nd_nof(f(:,:,:,iux+i-1),tmp(:,:,:,1),j) ! uij
call der_2nd_nof(f(:,:,:,iux+j-1),tmp(:,:,:,2),i) ! uji
!
! find off diagonal part of sij2 for full array
!
sij2=sij2+tmp(l1:l2,m1:m2,n1:n2,1)*tmp(l1:l2,m1:m2,n1:n2,2) &
+0.5*tmp(l1:l2,m1:m2,n1:n2,1)**2 &
+0.5*tmp(l1:l2,m1:m2,n1:n2,2)**2
enddo
enddo
!
! add diagonal part of sij2 for full array
!
tmp(:,:,:,1)=0
call der_2nd_nof(f(:,:,:,iux-1+1),tmp(:,:,:,1),1)
sij2=sij2+tmp(l1:l2,m1:m2,n1:n2,1)**2
call der_2nd_nof(f(:,:,:,iux-1+2),tmp(:,:,:,1),2)
sij2=sij2+tmp(l1:l2,m1:m2,n1:n2,1)**2
call der_2nd_nof(f(:,:,:,iux-1+3),tmp(:,:,:,1),3)
sij2=sij2+tmp(l1:l2,m1:m2,n1:n2,1)**2
!
f(l1:l2,m1:m2,n1:n2,ismagorinsky)=(C_smag*dxmax)**2.*sqrt(2*sij2)
!
! Diagnostics
!
if (ldiagnos) then
if (idiag_nu_LES /= 0) then
itype_name(idiag_nu_LES)=ilabel_sum
do mcount=m1,m2
do ncount=n1,n2
nu_smag=f(l1:l2,mcount,ncount,ismagorinsky)
call sum_mn_name(nu_smag,idiag_nu_LES)
enddo
enddo
endif
endif
!
!
!
!
endsubroutine calc_viscosity
!***********************************************************************
subroutine der_2nd_nof(var,tmp,j)
!
! 24-nov-03/nils: coded
!
use Cdata
!
real, dimension (mx,my,mz) :: var
real, dimension (mx,my,mz) :: tmp
integer :: j
!
intent (in) :: var,j
intent (out) :: tmp
!
tmp=0.
!
if (j==1 .and. nxgrid/=1) then
tmp( 1,:,:) = (- 3.*var(1,:,:) &
+ 4.*var(2,:,:) &
- 1.*var(3,:,:))/(2.*dx)
tmp(2:mx-1,:,:) = (- 1.*var(1:mx-2,:,:) &
+ 1.*var(3:mx ,:,:))/(2.*dx)
tmp( mx,:,:) = (+ 1.*var(mx-2,:,:) &
- 4.*var(mx-1,:,:) &
+ 3.*var(mx ,:,:))/(2.*dx)
endif
!
if (j==2 .and. nygrid/=1) then
tmp(:, 1,:) = (- 3.*var(:,1,:) &
+ 4.*var(:,2,:) &
- 1.*var(:,3,:))/(2.*dy)
tmp(:,2:my-1,:) = (- 1.*var(:,1:my-2,:) &
+ 1.*var(:,3:my ,:))/(2.*dy)
tmp(:, my,:) = (+ 1.*var(:,my-2,:) &
- 4.*var(:,my-1,:) &
+ 3.*var(:,my ,:))/(2.*dy)
endif
!
if (j==3 .and. nzgrid/=1) then
tmp(:,:, 1) = (- 3.*var(:,:,1) &
+ 4.*var(:,:,2) &
- 1.*var(:,:,3))/(2.*dz)
tmp(:,:,2:mz-1) = (- 1.*var(:,:,1:mz-2) &
+ 1.*var(:,:,3:mz ))/(2.*dz)
tmp(:,:, mz) = (+ 1.*var(:,:,mz-2) &
- 4.*var(:,:,mz-1) &
+ 3.*var(:,:,mz ))/(2.*dz)
endif
!
endsubroutine der_2nd_nof
!
!!**********************************************************************
subroutine calc_viscous_heat(f,df,glnrho,divu,rho1,cs2,TT1,shock,Hmax)
!
! calculate viscous heating term for right hand side of entropy equation
!
! 20-nov-02/tony: coded
!
use Cdata
use Mpicomm
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real, dimension (nx) :: rho1,TT1,cs2,Hmax
real, dimension (nx) :: sij2, divu,shock
real, dimension (nx,3) :: glnrho
!
! traceless strain matrix squared
!
call multm2_mn(sij,sij2)
!
select case (ivisc)
case default
if (headtt) print*,'no heating: ivisc=',ivisc
endselect
!
call keep_compiler_quiet(f)
call keep_compiler_quiet(cs2)
call keep_compiler_quiet(divu)
call keep_compiler_quiet(glnrho)
call keep_compiler_quiet(shock)
call keep_compiler_quiet(Hmax)
!
endsubroutine calc_viscous_heat
!***********************************************************************
subroutine calc_viscous_force(f,df,glnrho,divu,rho,rho1,shock,gshock,bij)
!
! calculate viscous heating term for right hand side of entropy equation
!
! 16-jul-04/nils: coded
!
use Cdata
use Mpicomm
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
real, dimension (mx,my,mz,mvar) :: df
real, dimension (nx,3,3) :: bij
real, dimension (nx,3) :: glnrho,del2u,del6u,graddivu,fvisc,sglnrho
real, dimension (nx,3) :: gshock,sgradnu_smag
real, dimension (nx,3) :: nusglnrho,tmp1,tmp2,gradnu_smag
real, dimension (nx) :: murho1,rho,rho1,divu,shock,SS12,sij2
integer :: i
!
intent (in) :: f, glnrho, rho1,rho
intent (out) :: df
!
! viscosity operator
! rho1 is pre-calculated in equ
!
select case (ivisc)
!
case ('smagorinsky')
!
! viscous force: nu_smag*(del2u+graddivu/3+2S.glnrho)+2S.gradnu_smag
! where nu_smag=(C_smag*dxmax)**2*sqrt(SS)
!
if (headtt) print*,'viscous force: Smagorinsky'
if (ldensity) then
call del2v_etc(f,iuu,del2u,GRADDIV=graddivu)
call multmv_mn(sij,glnrho,sglnrho)
call multsv_mn(f(l1:l2,m,n,ismagorinsky),sglnrho,nusglnrho)
tmp1=del2u+1./3.*graddivu
call multsv_mn(f(l1:l2,m,n,ismagorinsky),tmp1,tmp2)
call grad(f,ismagorinsky,gradnu_smag)
call multmv_mn(sij,gradnu_smag,sgradnu_smag)
fvisc=2*nusglnrho+tmp2+2*sgradnu_smag
diffus_nu=diffus_nu+f(l1:l2,m,n,ismagorinsky)*dxyz_2
endif
case ('smagorinsky_simplified')
if (headtt) print*, 'for ivisc=smagorinsky_simplified use visc_const'
!
! viscous force: nu_smag*(del2u+graddivu/3+2S.glnrho)
! where nu_smag=(C_smag*dxmax)**2*sqrt(SS)
!
if (headtt) print*,'viscous force: Smagorinsky_simplified'
if (ldensity) then
call del2v_etc(f,iuu,del2u,GRADDIV=graddivu)
call multmv_mn(sij,glnrho,sglnrho)
call multsv_mn(f(l1:l2,m,n,ismagorinsky),sglnrho,nusglnrho)
tmp1=del2u+1./3.*graddivu
call multsv_mn(f(l1:l2,m,n,ismagorinsky),tmp1,tmp2)
call grad(f,ismagorinsky,gradnu_smag)
call multmv_mn(2*sij,gradnu_smag,sgradnu_smag)
fvisc=2*nusglnrho+tmp2+sgradnu_smag
diffus_nu=diffus_nu+f(l1:l2,m,n,ismagorinsky)*dxyz_2
endif
case default
!
! Catch unknown values
!
if (lroot) print*, 'No such value for ivisc: ', trim(ivisc)
call stop_it('calc_viscous_forcing')
!
endselect
!
! Add ordinary viscosity if nu /= 0
!
if (nu /= 0.) then
!
! viscous force: nu*(del2u+graddivu/3+2S.glnrho)
! -- the correct expression for nu=const
!
fvisc=fvisc+2*nu*sglnrho+nu*(del2u+1./3.*graddivu)
diffus_nu=diffus_nu+nu*dxyz_2
endif
!
df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+fvisc
!
! set viscous time step
!
if (ldiagnos.and.idiag_dtnu/=0) then
call max_mn_name(diffus_nu/cdtv,idiag_dtnu,l_dt=.true.)
endif
!
call keep_compiler_quiet(divu)
call keep_compiler_quiet(shock)
call keep_compiler_quiet(gshock)
!
endsubroutine calc_viscous_force
!***********************************************************************
endmodule Viscosity