! $Id$
!
! This modules implements viscous heating and diffusion terms
! here for cases 1) nu constant, 2) mu = rho.nu 3) constant and
!
!** 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 0
!
! PENCILS PROVIDED fvisc(3); diffus_total; diffus_total2; diffus_total3
! PENCILS PROVIDED visc_heat; nu; gradnu(3), nu_smag, gnu_smag(3)
!
!***************************************************************
module Viscosity
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages
!
implicit none
!
include 'viscosity.h'
!
integer, parameter :: nvisc_max=4
character (len=labellen), dimension(nvisc_max) :: ivisc=''
character (len=labellen) :: lambda_profile='uniform'
real :: nu=0.0, nu_cspeed=0.5
real :: nu_tdep=0.0, nu_tdep_exponent=0.0, nu_tdep_t0=0.0, nu_tdep_toffset=0.0
real :: zeta=0.0, nu_mol=0.0, nu_hyper2=0.0, nu_hyper3=0.0
real :: nu_hyper3_mesh=5.0, nu_shock=0.0,nu_spitzer=0.0, nu_spitzer_max=0.0
real :: nu_jump=1.0, xnu=1.0, xnu2=1.0, znu=1.0, widthnu=0.1, widthnu2=0.1
real :: C_smag=0.0, gamma_smag=0.0, nu_jump2=1.0
real :: znu_shock=1.0, widthnu_shock=0.1, nu_jump_shock=1.0
real :: xnu_shock=1.0, dynu=0.0
real :: pnlaw=0.0, Lambda_V0=0.,Lambda_V1=0.,Lambda_H1=0.
real :: Lambda_V0t=0.,Lambda_V1t=0.,Lambda_V0b=0.,Lambda_V1b=0.
real :: rzero_lambda=impossible,wlambda=0.,rmax_lambda=impossible
real :: offamp_lambda=1.,r1_lambda=impossible,r2_lambda=impossible
real :: lambda_jump=0.,roffset_lambda=0.
real :: PrM_turb=0.0
real :: meanfield_nuB=0.0
real :: nu_infinity=0.,nu0=0.,non_newton_lambda=0.,carreau_exponent=0.
real :: nu_smag_Ma2_power=0.0
character (len=labellen) :: nnewton_type='none'
character (len=labellen) :: div_sld_visc='2nd'
real :: nnewton_tscale=0.0,nnewton_step_width=0.0
real, dimension(nx) :: xmask_vis=0, pnu=0.0
real, dimension(2) :: vis_xaver_range=(/-max_real,max_real/)
real, dimension(:), pointer :: etat_x, detat_x
real, dimension(:), pointer :: etat_y, detat_y
real, dimension(:), pointer :: etat_z, detat_z
real, dimension(3) :: nu_aniso_hyper3=0.0
real, dimension(mx) :: LV0_rprof,LV1_rprof,LH1_rprof,der_LV0_rprof,der_LV1_rprof
logical :: lvisc_first=.false.
logical :: lvisc_simplified=.false.
logical :: lvisc_nu_non_newtonian=.false.
logical :: lvisc_rho_nu_const=.false.
logical :: lvisc_rho_nu_const_bulk=.false.
logical :: lvisc_rho_nu_const_prefact=.false.
logical :: lvisc_sqrtrho_nu_const=.false.
logical :: lvisc_nu_cspeed=.false.
logical :: lvisc_mu_cspeed=.false.
logical :: lvisc_nu_const=.false.
logical :: lvisc_nu_tdep=.false.
logical :: lvisc_nu_tdep_t0_norm=.false.
logical :: lvisc_nu_prof=.false.
logical :: lvisc_nu_profx=.false.
logical :: lvisc_nu_profr=.false.
logical :: lvisc_nu_profr_powerlaw=.false.
logical :: lvisc_nu_profr_twosteps=.false.
logical :: lvisc_nu_profy_bound=.false.
logical :: lvisc_nut_from_magnetic=.false.
logical :: lvisc_nu_shock=.false.
logical :: lvisc_nu_shock_profz=.false.
logical :: lvisc_nu_shock_profr=.false.
logical :: lvisc_shock_simple=.false.
logical :: lvisc_hyper2_simplified=.false.
logical :: lvisc_hyper3_simplified=.false.
logical :: lvisc_hyper2_simplified_tdep=.false.
logical :: lvisc_hyper3_simplified_tdep=.false.
logical :: lvisc_hyper3_polar=.false.
logical :: lvisc_hyper3_mesh=.false.
logical :: lvisc_hyper3_mesh_residual=.false.
logical :: lvisc_hyper3_csmesh=.false.
logical :: lvisc_hyper3_rho_nu_const=.false.
logical :: lvisc_hyper3_mu_const_strict=.false.
logical :: lvisc_hyper3_nu_const_strict=.false.
logical :: lvisc_hyper3_cmu_const_strt_otf=.false.
logical :: lvisc_hyper3_rho_nu_const_symm=.false.
logical :: lvisc_hyper3_rho_nu_const_aniso=.false.
logical :: lvisc_hyper3_nu_const_aniso=.false.
logical :: lvisc_hyper3_rho_nu_const_bulk=.false.
logical :: lvisc_hyper3_nu_const=.false.
logical :: lvisc_smag_simplified=.false.
logical :: lvisc_smag_cross_simplified=.false.
logical :: lnusmag_as_aux=.false.
logical :: lvisc_snr_damp=.false.
logical :: lvisc_heat_as_aux=.false.
logical :: lvisc_forc_as_aux=.false.
logical :: lvisc_mixture=.false.
logical :: lvisc_spitzer=.false.
logical :: lvisc_slope_limited=.false.
logical :: limplicit_viscosity=.false.
logical :: lmeanfield_nu=.false.
logical :: lmagfield_nu=.false.
logical :: llambda_effect=.false.
logical :: luse_nu_rmn_prof=.false.
logical :: lvisc_smag_Ma=.false.
logical, pointer:: lviscosity_heat
logical :: lKit_Olem
logical :: lsld_notensor=.false.
logical :: lno_visc_heat_zbound=.false.
logical, pointer :: lcalc_uuavg
real :: no_visc_heat_z0=max_real, no_visc_heat_zwidth=0.0
real :: damp_sound=0.
real :: h_sld_visc=2.0, nlf_sld_visc=1.0
!
namelist /viscosity_run_pars/ &
limplicit_viscosity, nu, nu_tdep_exponent, &
lvisc_nu_tdep_t0_norm, nu_tdep_t0, nu_tdep_toffset, &
zeta, nu_hyper2, nu_hyper3, ivisc, nu_mol, C_smag, gamma_smag, nu_shock, &
nu_aniso_hyper3, lvisc_heat_as_aux,nu_jump,znu,xnu,xnu2,widthnu,widthnu2, &
pnlaw,llambda_effect,Lambda_V0,Lambda_V1,Lambda_H1, nu_hyper3_mesh, &
lambda_profile,rzero_lambda,wlambda,r1_lambda,r2_lambda,rmax_lambda,&
offamp_lambda,lambda_jump,lmeanfield_nu,lmagfield_nu,meanfield_nuB, &
PrM_turb, roffset_lambda, nu_spitzer, nu_jump2, nu_spitzer_max,&
widthnu_shock, znu_shock, xnu_shock, nu_jump_shock, dynu, &
nnewton_type,nu_infinity,nu0,non_newton_lambda,carreau_exponent,&
nnewton_tscale,nnewton_step_width,lKit_Olem,damp_sound,luse_nu_rmn_prof, &
h_sld_visc,nlf_sld_visc, lnusmag_as_aux, lsld_notensor, &
lvisc_smag_Ma, nu_smag_Ma2_power, nu_cspeed, lno_visc_heat_zbound, &
no_visc_heat_z0,no_visc_heat_zwidth, div_sld_visc ,lvisc_forc_as_aux, &
lvisc_rho_nu_const_prefact
!
! other variables (needs to be consistent with reset list below)
integer :: idiag_nu_tdep=0 ! DIAG_DOC: time-dependent viscosity
integer :: idiag_fviscm=0 ! DIAG_DOC: Mean value of viscous acceleration
integer :: idiag_fviscmin=0 ! DIAG_DOC: Min value of viscous acceleration
integer :: idiag_fviscmax=0 ! DIAG_DOC: Max value of viscous acceleration
integer :: idiag_fviscrmsx=0 ! DIAG_DOC: Rms value of viscous acceleration
! DIAG_DOC: for the vis_xaver_range
integer :: idiag_num=0 ! DIAG_DOC: Mean value of viscosity
integer :: idiag_nusmagm=0 ! DIAG_DOC: Mean value of Smagorinsky viscosity
integer :: idiag_nusmagmin=0 ! DIAG_DOC: Min value of Smagorinsky viscosity
integer :: idiag_nusmagmax=0 ! DIAG_DOC: Max value of Smagorinsky viscosity
integer :: idiag_nu_LES=0 ! DIAG_DOC: Mean value of Smagorinsky viscosity
integer :: idiag_visc_heatm=0 ! DIAG_DOC: Mean value of viscous heating
integer :: idiag_qfviscm=0 ! DIAG_DOC: $\left<\qv\cdot
! DIAG_DOC: \fv_{\rm visc}\right>$
integer :: idiag_ufviscm=0 ! DIAG_DOC: $\left<\uv\cdot
! DIAG_DOC: \fv_{\rm visc}\right>$
integer :: idiag_Sij2m=0 ! DIAG_DOC: $\left<\Strain^2\right>$
integer :: idiag_epsK=0 ! DIAG_DOC: $\left<2\nu\varrho\Strain^2\right>$
integer :: idiag_epsKint=0 ! DIAG_DOC: $\int(2\nu\varrho\Strain^2)\,dV$
integer :: idiag_epsK_LES=0 ! DIAG_DOC:
integer :: idiag_sijoiojm=0 ! DIAG_DOC: $\left<S_{i,j} \omega_i \omega_j\right>$
integer :: idiag_dtnu=0 ! DIAG_DOC: $\delta t/[c_{\delta t,{\rm v}}\,
! DIAG_DOC: \delta x^2/\nu_{\rm max}]$
! DIAG_DOC: \quad(time step relative to
! DIAG_DOC: viscous time step;
! DIAG_DOC: see \S~\ref{time-step})
integer :: idiag_dtnu3=0
integer :: idiag_meshRemax=0 ! DIAG_DOC: Max mesh Reynolds number
integer :: idiag_Reshock=0 ! DIAG_DOC: Mesh Reynolds number at shock
integer :: idiag_nuD2uxbxm=0 ! DIAG_DOC:
integer :: idiag_nuD2uxbym=0 ! DIAG_DOC:
integer :: idiag_nuD2uxbzm=0 ! DIAG_DOC:
!
! xy averaged diagnostics given in xyaver.in written every it1d timestep
!
integer :: idiag_fviscmz=0 ! XYAVG_DOC: $\left<2\nu\varrho u_i
! XYAVG_DOC: \mathcal{S}_{iz} \right>_{xy}$
! XYAVG_DOC: ($z$-component of viscous flux)
integer :: idiag_fviscsmmz=0 ! XYAVG_DOC: $\left<2\nu_{\rm Smag}\varrho u_i
! XYAVG_DOC: \mathcal{S}_{iz} \right>_{xy}$
! XYAVG_DOC: ($z$-component of viscous flux)
integer :: idiag_epsKmz=0 ! XYAVG_DOC: $\left<2\nu\varrho\Strain^2
! XYAVG_DOC: \right>_{xy}$
integer :: idiag_viscforcezmz=0 ! XYAVG_DOC: $\left<(\varrho\fv_{\rm visc})_z\right>_{xy}$
integer :: idiag_viscforcezupmz=0 ! XYAVG_DOC: $\left<(\varrho\fv_{\rm visc})_z\right>_{xy+}$
integer :: idiag_viscforcezdownmz=0 ! XYAVG_DOC: $\left<(\varrho\fv_{\rm visc})_z\right>_{xy-}$
!
! yz averaged diagnostics given in yzaver.in written every it1d timestep
!
integer :: idiag_fviscmx=0 ! YZAVG_DOC: $\left<2\nu\varrho u_i
! YZAVG_DOC: \mathcal{S}_{ix} \right>_{yz}$
! YZAVG_DOC: ($x$-component of viscous flux)
integer :: idiag_numx=0 ! YZAVG_DOC: $\left< \nu \right>_{yz}$
! YZAVG_DOC: ($yz$-averaged viscosity)
!
! z averaged diagnostics given in zaver.in
!
integer :: idiag_fviscmxy=0 ! ZAVG_DOC: $\left<2\nu\varrho u_i
! ZAVG_DOC: \mathcal{S}_{ix} \right>_{z}$
! ZAVG_DOC: ($x$-xomponent of viscous flux)
integer :: idiag_fviscsmmxy=0 ! ZAVG_DOC: $\left<2\nu_{\rm Smag}\varrho u_i
! ZAVG_DOC: \mathcal{S}_{ix} \right>_{z}$
! ZAVG_DOC: ($x$-xomponent of viscous flux)
integer :: idiag_fviscymxy=0 ! ZAVG_DOC: $\left<2\nu\varrho u_i
! ZAVG_DOC: \mathcal{S}_{iy} \right>_{z}$
! ZAVG_DOC: ($y$-xomponent of viscous flux)
!
! phi averaged diagnostics given in phiaver.in
!
integer :: idiag_fviscrsphmphi ! PHIAVG_DOC: $\left<2\nu\varrho u_i
! PHIAVG_DOC: \mathcal{S}_{ir} \right>_\varphi$
! PHIAVG_DOC: ($r$-xomponent of viscous flux)
!
! Module Variables
!
real, dimension(mz) :: eth0z = 0.0
!
contains
!***********************************************************************
subroutine register_viscosity
use FArrayManager, only: farray_register_auxiliary
!
! 19-nov-02/tony: coded
! 30-sep-15/Joern+MR: changes for slope-limited diffusion
! 03-apr-20/joern: restructured and fixed slope-limited diffusion
!
! Identify version number.
!
if (lroot) call svn_id( &
"$Id$")
!
! Register characteristic speed: sld_char as auxilliary variable
! Needed for slope limited diffusion
!
if (any(ivisc=='nu-slope-limited')) then
lslope_limit_diff = .true.
if (dimensionality<3) lisotropic_advection=.true.
if (isld_char == 0) then
call farray_register_auxiliary('sld_char',isld_char,communicated=.true.)
if (lroot) write(15,*) 'sld_char = fltarr(mx,my,mz)*one'
aux_var(aux_count)=',sld_char'
if (naux+naux_com < maux+maux_com) aux_var(aux_count)=trim(aux_var(aux_count))//' $'
aux_count=aux_count+1
endif
endif
!
! Register nusmag as auxilliary variable
!
if (lnusmag_as_aux.and.any(ivisc=='smagorinsky')) then
call farray_register_auxiliary('nusmag',inusmag,communicated=.true.)
if (lroot) write(15,*) 'nusmag = fltarr(mx,my,mz)*one'
aux_var(aux_count)=',nusmag'
if (naux+naux_com < maux+maux_com) aux_var(aux_count)=trim(aux_var(aux_count))//' $'
aux_count=aux_count+1
endif
!
! Register an extra aux slot for dissipation rate if requested (so
! visc_heat is written to snapshots and can be easily analyzed later).
!
if (lvisc_heat_as_aux) then
call farray_register_auxiliary('visc_heat',ivisc_heat)
if (lroot) write(15,*) 'visc_heat = fltarr(mx,my,mz)*one'
aux_var(aux_count)=',visc_heat'
if (naux+naux_com < maux+maux_com) aux_var(aux_count)=trim(aux_var(aux_count))//' $'
aux_count=aux_count+1
endif
!
! Register an 3 extra aux slot for viscouse force (accelaration) if requested (so
! visc_forc is written to snapshots and can be easily analyzed later).
!
if (lvisc_forc_as_aux) then
call farray_register_auxiliary('visc_forc',ivisc_forc,vector=3)
if (lroot) write(15,*) 'visc_forc = fltarr(mx,my,mz,3)*one'
aux_var(aux_count)=',visc_forc'
if (naux+naux_com < maux+maux_com) aux_var(aux_count)=trim(aux_var(aux_count))//' $'
aux_count=aux_count+3
ivisc_forcx=ivisc_forc;ivisc_forcy=ivisc_forc+1;ivisc_forcz=ivisc_forc+2
endif
!
endsubroutine register_viscosity
!***********************************************************************
subroutine initialize_viscosity
!
! 20-nov-02/tony: coded
!
use EquationOfState, only: get_stratz
use Mpicomm, only: stop_it
use SharedVariables, only: put_shared_variable,get_shared_variable
use Sub, only: write_zprof, write_yprof, step
!
integer :: i
integer :: ierr
!
! Default viscosity.
!
if ( (nu/=0.0).and.(ivisc(1)=='') ) ivisc(1)='nu-const'
!
! Some viscosity types need the rate-of-strain tensor and grad(lnrho)
!
lvisc_simplified=.false.
lvisc_rho_nu_const=.false.
lvisc_rho_nu_const_bulk=.false.
lvisc_sqrtrho_nu_const=.false.
lvisc_nu_cspeed=.false.
lvisc_mu_cspeed=.false.
lvisc_nu_const=.false.
lvisc_nu_tdep=.false.
lvisc_nu_prof=.false.
lvisc_nu_profx=.false.
lvisc_nu_profr=.false.
lvisc_nu_profr_powerlaw=.false.
lvisc_nu_profr_twosteps=.false.
lvisc_nu_profy_bound=.false.
lvisc_nut_from_magnetic=.false.
lvisc_nu_shock=.false.
lvisc_nu_shock_profz=.false.
lvisc_nu_shock_profr=.false.
lvisc_shock_simple=.false.
lvisc_hyper2_simplified=.false.
lvisc_hyper3_simplified=.false.
lvisc_hyper2_simplified_tdep=.false.
lvisc_hyper3_simplified_tdep=.false.
lvisc_hyper3_polar=.false.
lvisc_hyper3_mesh=.false.
lvisc_hyper3_mesh_residual=.false.
lvisc_hyper3_csmesh=.false.
lvisc_hyper3_rho_nu_const=.false.
lvisc_hyper3_rho_nu_const_symm=.false.
lvisc_hyper3_mu_const_strict=.false.
lvisc_hyper3_nu_const_strict=.false.
lvisc_hyper3_cmu_const_strt_otf=.false.
lvisc_hyper3_rho_nu_const_aniso=.false.
lvisc_hyper3_nu_const_aniso=.false.
lvisc_hyper3_rho_nu_const_bulk=.false.
lvisc_hyper3_nu_const=.false.
lvisc_smag=.false.
lvisc_smag_simplified=.false.
lvisc_smag_cross_simplified=.false.
lvisc_snr_damp=.false.
lvisc_spitzer=.false.
lvisc_slope_limited=.false.
!
do i=1,nvisc_max
select case (ivisc(i))
case ('nu-simplified','simplified', '0')
if (lroot) print*,'viscous force: nu*del2v'
lvisc_simplified=.true.
case ('nu-non-newtonian')
if (lroot) print*,'viscous force: div(nu*Sij)'
lvisc_nu_non_newtonian=.true.
case ('rho-nu-const','rho_nu-const', '1')
if (lroot) print*,'viscous force: mu/rho*(del2u+graddivu/3)'
lvisc_rho_nu_const=.true.
case ('rho-nu-const-bulk')
if (lroot) print*,'viscous force: (zeta/rho)*graddivu'
lvisc_rho_nu_const_bulk=.true.
case ('sqrtrho-nu-const')
if (lroot) print*,'viscous force: mu/sqrt(rho)*(del2u+graddivu/3+S.glnrho)'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_sqrtrho_nu_const=.true.
case ('nu-cspeed')
if (lroot) print*,'viscous force: mu*sqrt(TT)*(del2u+graddivu/3+2S.glnrho)'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_cspeed=.true.
case ('mu-cspeed')
if (lroot) print*,'viscous force: nu*sqrt(TT)/rho*(del2u+graddivu/3)'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_mu_cspeed=.true.
case ('nu-const')
if (lroot) print*,'viscous force: nu*(del2u+graddivu/3+2S.glnrho)'
if (nu/=0.) lpenc_requested(i_sij)=.true.
! if (meanfield_nuB/=0.) lpenc_requested(i_b2)=.true.
lpenc_requested(i_glnrho)=.true.
lvisc_nu_const=.true.
case ('nu-tdep')
if (lroot) print*,'time-dependent nu*(del2u+graddivu/3+2S.glnrho)'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_tdep=.true.
case ('nu-prof')
if (lroot) print*,'viscous force with a vertical profile for nu'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_prof=.true.
case ('nu-profx')
if (lroot) print*,'viscous force with a horizontal profile for nu'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_profx=.true.
case ('nu-profr')
if (lroot) print*,'viscous force with a radial profile for nu'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_profr=.true.
case ('nu-profr-twosteps')
if (lroot) print*,'viscous force with a radial profile for nu'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_profr_twosteps=.true.
case ('nu-profr-powerlaw','powerlaw','power-law')
if (lroot) print*,'viscous force with a power law profile'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_profr_powerlaw=.true.
case ('nu-profy-bound')
if (lroot) print*,'viscous force with a horizontal profile for nu'
if (lroot) print*,'where the nu is enhanced near the y boundaries'
if (nu/=0.) lpenc_requested(i_sij)=.true.
lvisc_nu_profy_bound=.true.
case ('nut-from-magnetic')
if (lroot) print*,'nut-from-magnetic via shared variables'
if (PrM_turb/=0.) lpenc_requested(i_sij)=.true.
lvisc_nut_from_magnetic=.true.
case ('nu-shock','shock')
if (lroot) print*,'viscous force: nu_shock*(XXXXXXXXXXX)'
lvisc_nu_shock=.true.
if (.not. lshock) &
call stop_it('initialize_viscosity: shock viscosity'// &
' but module setting SHOCK=noshock')
case ('nu-shock-profz')
if (lroot) print*,'viscous force: nu_shock*(XXXXXXXXXXX) with a vertical profile'
lvisc_nu_shock_profz=.true.
if (.not. lshock) &
call stop_it('initialize_viscosity: shock viscosity'// &
' but module setting SHOCK=noshock')
case ('nu-shock-profr')
if (lroot) print*,'viscous force: nu_shock*(XXXXXXXXXXX) with a radial profile'
lvisc_nu_shock_profr=.true.
if (.not. lshock) &
call stop_it('initialize_viscosity: shock viscosity'// &
' but module setting SHOCK=noshock')
case ('shock_simple', 'shock-simple')
if (lroot) print *, 'viscous force: div(nu_shock*grad(uu_i)))'
lvisc_shock_simple = .true.
if (.not. lshock) call stop_it('initialize_viscosity: a SHOCK module is required. ')
case ('hyper2-simplified','hyper2_simplified', 'hyper4')
if (lroot) print*,'viscous force: -nu_hyper*del4v'
lvisc_hyper2_simplified=.true.
case ('hyper3-simplified','hyper3_simplified', 'hyper6')
if (lroot) print*,'viscous force: nu_hyper*del6v'
lvisc_hyper3_simplified=.true.
case ('hyper2-simplified-tdep')
if (lroot) print*,'viscous force: -nu*del4v'
lvisc_hyper2_simplified_tdep=.true.
case ('hyper3-simplified-tdep')
if (lroot) print*,'viscous force: nu*del6v'
lvisc_hyper3_simplified_tdep=.true.
case ('hyper3-cyl','hyper3_cyl','hyper3-sph','hyper3_sph')
if (lroot) print*,'viscous force: nu_hyper3/pi^4 *(Deltav)^6/Deltaq^2'
lvisc_hyper3_polar=.true.
case ('hyper3-mesh','hyper3_mesh')
if (lroot) print*,'viscous force: nu_hyper3_mesh/pi^5 *(Deltav)^6/Deltaq'
lvisc_hyper3_mesh=.true.
case ('hyper3-mesh-residual','hyper3_mesh-residual','hyper3-mesh_residual','hyper3_mesh_residual')
if (lroot) print*,'viscous force: nu_hyper3_mesh/pi^5 *(Delta(v-<v>)^6/Deltaq'
lvisc_hyper3_mesh_residual=.true.
call get_shared_variable('lcalc_uuavg',lcalc_uuavg,ierr)
if (ierr/=0) call fatal_error("initialize_viscosity","could not get lcalc_uuavg")
lcalc_uuavg=.true.
case ('hyper3-csmesh')
if (lroot) print*,'viscous force: c_s*nu_hyper3_mesh/pi^5 *(Deltav)^6/Deltaq'
lvisc_hyper3_csmesh=.true.
case ('hyper3-rho-nu-const','hyper3_rho_nu-const')
if (lroot) print*,'viscous force: nu_hyper/rho*del6v'
lvisc_hyper3_rho_nu_const=.true.
case ('hyper3-rho-nu-const-symm','hyper3_rho_nu-const_symm')
if (lroot) print*,'viscous force(i): nu_hyper/rho*(del6ui+der5(divu,i))'
lvisc_hyper3_rho_nu_const_symm=.true.
case ('hyper3-mu-const-strict','hyper3_mu-const_strict')
if (lroot) print*, 'viscous force(i): '// &
'nu_hyper/rho*(del2(del2(del2(u)))+del2(del2(grad(divu))))'
if (.not.lhyperviscosity_strict) &
call stop_it('initialize_viscosity: This viscosity type'//&
' cannot be used with HYPERVISC_STRICT=nohypervisc_strict')
lvisc_hyper3_mu_const_strict=.true.
case ('hyper3-mu-strict-onthefly')
if (lroot) print*, 'viscous force(i): '// &
'nu_hyper/rho*(del2(del2(del2(u)))+del2(del2(grad(divu))))'
lvisc_hyper3_cmu_const_strt_otf=.true.
case ('hyper3-nu-const-strict','hyper3_nu-const_strict')
if (lroot) print*, 'viscous force(i): 1/rho*div[2*rho*nu_3*S^(3)]'
if (.not.lhyperviscosity_strict) &
call stop_it('initialize_viscosity: This viscosity type'//&
' cannot be used with HYPERVISC_STRICT=nohypervisc_strict')
lvisc_hyper3_nu_const_strict=.true.
case ('hyper3-rho-nu-const-aniso','hyper3_rho_nu-const_aniso')
if (lroot) print*,&
'viscous force(i): 1/rho*(nu.del6)ui'
lvisc_hyper3_rho_nu_const_aniso=.true.
case ('hyper3-nu-const-aniso','hyper3_nu-const_aniso')
if (lroot) print*,&
'viscous force(i): (nu.del6)ui + ((nu.uij5).glnrho)'
lpenc_requested(i_uij5)=.true.
lpenc_requested(i_glnrho)=.true.
lvisc_hyper3_nu_const_aniso=.true.
case ('hyper3-rho-nu-const-bulk','hyper3_rho_nu-const_bulk')
if (lroot) print*,'viscous force: duj/dt = nu_hyper/rho*d6uj/dxj6'
lvisc_hyper3_rho_nu_const_bulk=.true.
case ('hyper3-nu-const','hyper3_nu-const','hyper3-nu_const')
if (lroot) print*,'viscous force: nu*(del6u+S.glnrho)'
lpenc_requested(i_uij5)=.true.
lvisc_hyper3_nu_const=.true.
case ('smagorinsky')
if (lroot) print*,'viscous force: Smagorinsky'
lpenc_requested(i_sij)=.true.
lvisc_smag=.true.
case ('smagorinsky-simplified','smagorinsky_simplified')
if (lroot) print*,'viscous force: Smagorinsky_simplified'
lpenc_requested(i_sij)=.true.
lvisc_smag_simplified=.true.
case ('smagorinsky-cross-simplif','smagorinsky_cross_simplif')
if (lroot) print*,'viscous force: Smagorinsky_simplified'
lvisc_smag_cross_simplified=.true.
! case ('snr-damp','snr_damp')
! if (lroot) print*,'viscous force: SNR damping'
! lvisc_snr_damp=.true.
case ('nu-mixture')
if (lroot) print*,'viscous force: nu is calculated for a mixture'
lvisc_mixture=.true.
case ('nu-spitzer')
if (lroot) print*,'viscous force: temperature dependent nu'
lpenc_requested(i_sij)=.true.
lvisc_spitzer=.true.
case ('nu-slope-limited')
if (lroot) print*,'viscous force: slope-limited diffusion'
if (lroot) print*,'viscous force: using ',trim(div_sld_visc),' order'
lvisc_slope_limited=.true.
case ('none',' ')
! do nothing
case default
if (lroot) print*, 'No such value for ivisc(',i,'): ', trim(ivisc(i))
call stop_it('calc_viscous_forcing')
endselect
enddo
!
! If we're timestepping, die or warn if the viscosity coefficient that
! corresponds to the chosen viscosity type is not set.
!
if (lrun) then
if (lroot) then
if ((lvisc_simplified.or.lvisc_rho_nu_const.or. &
lvisc_sqrtrho_nu_const.or.lvisc_nu_const.or. &
lvisc_nu_tdep.or.lvisc_nu_cspeed.or.&
lvisc_mu_cspeed).and.nu==0.0) &
call warning('initialize_viscosity', &
'Viscosity coefficient nu is zero!')
if ((lvisc_rho_nu_const_bulk).and.zeta==0.0) &
call warning('initialize_viscosity', &
'Viscosity coefficient zeta is zero!')
endif
if (lvisc_hyper2_simplified.and.nu_hyper2==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_hyper2 is zero!')
if ( (lvisc_hyper3_simplified.or.lvisc_hyper3_rho_nu_const.or. &
lvisc_hyper3_rho_nu_const_bulk.or.lvisc_hyper3_nu_const.or. &
lvisc_hyper3_rho_nu_const_symm.or. &
lvisc_hyper3_polar.or.&
lvisc_hyper3_mu_const_strict .or. &
lvisc_hyper3_nu_const_strict .or. &
lvisc_hyper3_cmu_const_strt_otf ).and. &
nu_hyper3==0.0 ) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_hyper3 is zero!')
if (lvisc_hyper3_mesh.and.nu_hyper3_mesh==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_hyper3_mesh is zero!')
if (lvisc_hyper3_mesh_residual.and.nu_hyper3_mesh==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_hyper3_mesh_residual is zero!')
if (lvisc_hyper3_csmesh.and.nu_hyper3_mesh==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_hyper3_mesh is zero!')
if (lvisc_spitzer.and.nu_spitzer==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_spitzer is zero!')
if ( (lvisc_hyper3_rho_nu_const_aniso.or.lvisc_hyper3_nu_const_aniso).and.&
((nu_aniso_hyper3(1)==0. .and. nxgrid/=1 ).or. &
(nu_aniso_hyper3(2)==0. .and. nygrid/=1 ).or. &
(nu_aniso_hyper3(3)==0. .and. nzgrid/=1 )) ) &
call fatal_error('initialize_viscosity', &
'A viscosity coefficient of nu_aniso_hyper3 is zero!')
if ( (lvisc_smag.or.lvisc_smag_simplified.or.lvisc_smag_cross_simplified).and. &
C_smag==0.0 ) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient C_smag is zero!')
if (lvisc_nu_shock .and.nu_shock==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_shock is zero!')
if (lvisc_nu_shock_profz .and.nu_shock==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_shock is zero!')
if (lvisc_nu_shock_profr .and.nu_shock==0.0) &
call fatal_error('initialize_viscosity', &
'Viscosity coefficient nu_shock is zero!')
if (lvisc_shock_simple .and. nu_shock == 0.0) call fatal_error('initialize_viscosity', 'nu_shock is zero. ')
!
! Dynamical hyper-diffusivity operates only for mesh formulation of hyper-viscosity
!
if (ldynamical_diffusion.and. &
.not.(lvisc_hyper3_mesh.or.lvisc_hyper3_mesh_residual.or.lvisc_hyper3_csmesh)) then
call fatal_error("initialize_viscosity",&
"Dynamical diffusion requires mesh hyper-diffusion, switch ivisc='hyper3-mesh' "// &
"'hyper3-mesh-residual', or 'hyper3-csmesh'")
endif
if (.not.ldensity) then
if (lvisc_smag.or.lvisc_smag_simplified.or.lvisc_smag_cross_simplified) &
call warning('initialize_viscosity', &
"ldensity better be .true. for ivisc='smagorinsky*'")
endif
endif
if (lyinyang) then
if (lvisc_nu_prof.or.lvisc_nu_shock_profz.or.lvisc_nut_from_magnetic) &
call fatal_error("initialize_viscosity",&
"z dependent profiles not implemented for Yin-Yang grid")
endif
if (lvisc_nu_prof) then
!
! Write out viscosity z-profile.
! At present only correct for Cartesian geometry
! The actually profile is generated below and stored in pnu.
!
call write_zprof('visc', nu + (nu*(nu_jump-1.))*step(z(n1:n2),znu,-widthnu))
endif
if (lvisc_nu_profy_bound) then
!
! Write out viscosity y-profile
! The actually profile is generated below and stored in pnu.
!
call write_yprof('visc_profy_bound', nu + (nu*(nu_jump-1.))* &
(step(y(m1:m2),xyz1(2)-3*dynu, dynu) + step(y(m1:m2),xyz0(2)+3*dynu, -dynu)))
endif
if (lvisc_nu_shock_profz .or. lvisc_nu_shock_profr) then
!
! Write out shock viscosity z-profile.
! At present only correct for Cartesian geometry
!
if (lvisc_nu_shock_profz) &
call write_zprof('visc_shock', &
nu_shock+(nu_shock*(nu_jump_shock-1.))*step(z(n1:n2),znu_shock,-widthnu_shock))
!
! Write out shock viscosity r-profile.
! At present only correct for spherical and cylindrical geometry.
!
if (lvisc_nu_shock_profr) &
call write_zprof('visc_shock', &
nu_shock+(nu_shock*(nu_jump_shock-1.))*step(x(l1:l2),xnu_shock,-widthnu_shock))
endif
!
lnusmag_as_aux = lnusmag_as_aux.and.lvisc_smag
!
! Shared variables.
!
call put_shared_variable('lvisc_hyper3_nu_const_strict',lvisc_hyper3_nu_const_strict, &
caller='initialize_viscosity')
call put_shared_variable('nu',nu)
call get_shared_variable('lviscosity_heat',lviscosity_heat)
call put_shared_variable('llambda_effect',llambda_effect)
!
! For Boussinesq, density is constant, so must use lvisc_simplified.
!
if (lrun.and.lboussinesq) then
if (.not.lvisc_simplified) call fatal_error("initialize_viscosity", &
"Boussinesq works only if ivisc='simplified' in run.in")
if (lroot) print*,'viscosity: viscous force nu*del2v with nu=',nu
endif
!
! DM Note:
! initialization of lambda should come after sharing the llambda_effect
! variable because even if there is no lambda_effect that has to be
! communicated to the boundary conditions too. In case llambda_effect is true
! more variables are shared inside initialize_lambda.
! initialize lambda effect.
!
if (llambda_effect) call initialize_lambda
!
! Check for possibility of getting etat profile and gradient from
! the magnetic meanfield module.
!
if (PrM_turb/=0.) then
if (lmagn_mf) then
call get_shared_variable('etat_x',etat_x)
call get_shared_variable('etat_y',etat_y)
call get_shared_variable('etat_z',etat_z)
call get_shared_variable('detat_x',detat_x)
call get_shared_variable('detat_y',detat_y)
call get_shared_variable('detat_z',detat_z)
do n=n1,n2
print*,ipz,z(n),etat_z(n),detat_z(n)
enddo
endif
endif
!
!
! Compute mask for x-averaging where x is in vis_xaver_range.
! Normalize such that the average over the full domain
! gives still unity.
!
if (l1 == l2) then
xmask_vis = 1.
else
where (x(l1:l2) >= vis_xaver_range(1) .and. x(l1:l2) <= vis_xaver_range(2))
xmask_vis = 1.
elsewhere
xmask_vis = 0.
endwhere
vis_xaver_range(1) = max(vis_xaver_range(1), xyz0(1))
vis_xaver_range(2) = min(vis_xaver_range(2), xyz1(1))
if (lspherical_coords) then
xmask_vis = xmask_vis * (xyz1(1)**3 - xyz0(1)**3) &
/ (vis_xaver_range(2)**3 - vis_xaver_range(1)**3)
elseif (lcylindrical_coords) then
xmask_vis = xmask_vis * (xyz1(1)**2 - xyz0(1)**2) &
/ (vis_xaver_range(2)**2 - vis_xaver_range(1)**2)
else
xmask_vis = xmask_vis*Lxyz(1) &
/ (vis_xaver_range(2) - vis_xaver_range(1))
endif
endif
!
! Get background energy stratification, if any.
!
if (lstratz .and. lthermal_energy) call get_stratz(z, eth0z=eth0z)
!
! debug output
!
if (lroot.and.ip<14) print*,'xmask_vis=',xmask_vis
!
endsubroutine initialize_viscosity
!***********************************************************************
subroutine initialize_lambda
!
! DM 26-05-2010: cut out of intialize viscosity
!
use Sub, only: step,der_step,stepdown,der_stepdown
use SharedVariables, only: put_shared_variable,get_shared_variable
!
if ((Lambda_V0==0).and.(Lambda_V1==0).and.(Lambda_H1==0)) &
call warning('initialize_lambda', &
'You have chose llambda_effect=T but, all Lambda coefficients to be zero!')
if ((Lambda_V0==0).and.((Lambda_V1/=0).or.(Lambda_H1==0))) &
call warning('initialize_lambda', &
'Lambda effect: V_zero=0 but V1 or H1 nonzero')
!
! Select the profile of Lambda, default is uniform. At present (May 2010) the
! only other coded profile is radial step.
!
select case (lambda_profile)
case ('radial_step_V0')
if (lroot) print*,'lambda profile radial_step_V0, rzero_lambda, wlambda:',&
rzero_lambda,wlambda
LV0_rprof=step(x,rzero_lambda,wlambda)
der_LV0_rprof=der_step(x,rzero_lambda,wlambda)
LV1_rprof=1.;LH1_rprof=1.
der_LV1_rprof=0.
case ('radial_step')
if (lroot) print*,'lambda profile radial_step, rzero_lambda, wlambda:',&
rzero_lambda,wlambda
LV0_rprof=step(x,rzero_lambda,wlambda)
LV1_rprof=step(x,rzero_lambda,wlambda)
LH1_rprof=step(x,rzero_lambda,wlambda)
der_LV0_rprof=der_step(x,rzero_lambda,wlambda)
der_LV1_rprof=der_step(x,rzero_lambda,wlambda)
case ('top_hat')
if (lroot) print*,'lambda profile top_hat, rzero_lambda, rmax_lambda,roffset_lambda, wlambda:',&
rzero_lambda,rmax_lambda,roffset_lambda,wlambda
LV0_rprof=step(x,rzero_lambda+roffset_lambda,wlambda) &
-step(x,rmax_lambda+roffset_lambda,wlambda)
LV1_rprof=step(x,rzero_lambda,wlambda)-step(x,rmax_lambda,wlambda)
LH1_rprof=step(x,rzero_lambda,wlambda)-step(x,rmax_lambda,wlambda)
der_LV0_rprof=der_step(x,rzero_lambda+roffset_lambda,wlambda) &
-der_step(x,rmax_lambda+roffset_lambda,wlambda)
der_LV1_rprof=der_step(x,rzero_lambda,wlambda) &
-der_step(x,rmax_lambda,wlambda)
case ('V1H1_roff')
LV0_rprof=1.;der_LV0_rprof=0.
LV1_rprof=1.+offamp_lambda*stepdown(x,rmax_lambda,wlambda)
LH1_rprof=1.+offamp_lambda*stepdown(x,rmax_lambda,wlambda)
der_LV1_rprof=offamp_lambda*der_stepdown(x,rmax_lambda,wlambda)
if (lroot) print*,'LV1_rprof',LV1_rprof
if (lroot) print*,'LH1_rprof',LH1_rprof
case ('roff')
LV0_rprof=1.+stepdown(x,rmax_lambda,wlambda)
der_LV0_rprof=1.+der_stepdown(x,rmax_lambda,wlambda)
LV1_rprof=1.+stepdown(x,rmax_lambda,wlambda)
LH1_rprof=1.+stepdown(x,rmax_lambda,wlambda)
der_LV1_rprof=der_stepdown(x,rmax_lambda,wlambda)
if (lroot) print*,'LV1_rprof',LV1_rprof
if (lroot) print*,'LH1_rprof',LH1_rprof
case ('uniform')
if (lroot) print*,'lambda profile :',lambda_profile
LV0_rprof=1.;LV1_rprof=1;LH1_rprof=1.
der_LV0_rprof=0.;der_LV1_rprof=0
case ('V0jump')
if (lroot) print*,'lambda_profile, radial_step, rzero_lambda, wlambda:',&
lambda_profile,rzero_lambda,wlambda
LV0_rprof=1.+(lambda_jump/Lambda_V0)*step(x,rzero_lambda,wlambda)
der_LV0_rprof=(lambda_jump/Lambda_V0)*der_step(x,rzero_lambda,wlambda)
LV1_rprof=1;LH1_rprof=1.
der_LV1_rprof=0
case default
call fatal_error('initialize_lambda',&
'default lambda_profile is uniform ! ')
endselect
lambda_V0t=lambda_V0*LV0_rprof(nx)
lambda_V1t=lambda_V1*LV1_rprof(nx)
lambda_V0b=lambda_V0*LV0_rprof(1)
lambda_V1b=lambda_V1*LV1_rprof(1)
call put_shared_variable('Lambda_V0t',Lambda_V0t,caller='initialize_lambda')
call put_shared_variable('Lambda_V1t',Lambda_V1t)
call put_shared_variable('Lambda_V0b',Lambda_V0b)
call put_shared_variable('Lambda_V1b',Lambda_V1b)
call put_shared_variable('Lambda_H1',Lambda_H1)
call put_shared_variable('LH1_rprof',LH1_rprof)
!
endsubroutine initialize_lambda
!***********************************************************************
subroutine read_viscosity_run_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=viscosity_run_pars, IOSTAT=iostat)
!
endsubroutine read_viscosity_run_pars
!***********************************************************************
subroutine write_viscosity_run_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=viscosity_run_pars)
!
endsubroutine write_viscosity_run_pars
!***********************************************************************
subroutine rprint_viscosity(lreset,lwrite)
!
! Writes ishock to index.pro file
!
! 24-nov-03/tony: adapted from rprint_ionization
!
use Diagnostics, only: parse_name
!
logical :: lreset
logical, intent(in), optional :: lwrite
integer :: iname,inamex,inamez,ixy,irz
!
! reset everything in case of reset
! (this needs to be consistent with what is defined above!)
!
if (lreset) then
idiag_dtnu=0; idiag_dtnu3=0; idiag_nu_LES=0; idiag_Sij2m=0
idiag_epsK=0; idiag_epsKint=0; idiag_epsK_LES=0; idiag_sijoiojm=0
idiag_visc_heatm=0; idiag_meshRemax=0; idiag_Reshock=0
idiag_nuD2uxbxm=0; idiag_nuD2uxbym=0; idiag_nuD2uxbzm=0
idiag_nu_tdep=0; idiag_fviscm=0 ; idiag_fviscrmsx=0
idiag_fviscmz=0; idiag_fviscmx=0; idiag_fviscmxy=0
idiag_epsKmz=0; idiag_numx=0; idiag_fviscymxy=0
idiag_fviscsmmz=0; idiag_fviscsmmxy=0; idiag_ufviscm=0
idiag_fviscmax=0; idiag_fviscmin=0; idiag_fviscrsphmphi=0
idiag_viscforcezmz=0; idiag_viscforcezupmz=0; idiag_viscforcezdownmz=0
endif
!
! iname runs through all possible names that may be listed in print.in
!
if (lroot.and.ip<14) print*,'rprint_viscosity: run through parse list'
do iname=1,nname
call parse_name(iname,cname(iname),cform(iname),'nu_tdep',idiag_nu_tdep)
call parse_name(iname,cname(iname),cform(iname),'fviscm',idiag_fviscm)
call parse_name(iname,cname(iname),cform(iname),'fviscmin',idiag_fviscmin)
call parse_name(iname,cname(iname),cform(iname),'qfviscm',idiag_qfviscm)
call parse_name(iname,cname(iname),cform(iname),'ufviscm',idiag_ufviscm)
call parse_name(iname,cname(iname),cform(iname),'fviscmax',idiag_fviscmax)
call parse_name(iname,cname(iname),cform(iname),'fviscrmsx',idiag_fviscrmsx)
call parse_name(iname,cname(iname),cform(iname),'num',idiag_num)
call parse_name(iname,cname(iname),cform(iname),'nusmagm',idiag_nusmagm)
call parse_name(iname,cname(iname),cform(iname),'nusmagmin',idiag_nusmagmin)
call parse_name(iname,cname(iname),cform(iname),'nusmagmax',idiag_nusmagmax)
call parse_name(iname,cname(iname),cform(iname),'dtnu',idiag_dtnu)
call parse_name(iname,cname(iname),cform(iname),'dtnu3',idiag_dtnu3)
call parse_name(iname,cname(iname),cform(iname),'nu_LES',idiag_nu_LES)
call parse_name(iname,cname(iname),cform(iname),'visc_heatm',idiag_visc_heatm)
call parse_name(iname,cname(iname),cform(iname),'Sij2m',idiag_Sij2m)
call parse_name(iname,cname(iname),cform(iname),'sijoiojm',idiag_sijoiojm)
call parse_name(iname,cname(iname),cform(iname),'epsK',idiag_epsK)
call parse_name(iname,cname(iname),cform(iname),'epsKint',idiag_epsKint)
call parse_name(iname,cname(iname),cform(iname),'epsK_LES',idiag_epsK_LES)
call parse_name(iname,cname(iname),cform(iname),'meshRemax',idiag_meshRemax)
call parse_name(iname,cname(iname),cform(iname),'Reshock',idiag_Reshock)
enddo
!
! Check for those quantities for which we want xy-averages.
!
do inamez=1,nnamez
call parse_name(inamez,cnamez(inamez),cformz(inamez),'fviscmz',idiag_fviscmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'fviscsmmz',idiag_fviscsmmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'epsKmz',idiag_epsKmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'viscforcezmz',idiag_viscforcezmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'viscforcezupmz',idiag_viscforcezupmz)
call parse_name(inamez,cnamez(inamez),cformz(inamez),'viscforcezdownmz',idiag_viscforcezdownmz)
enddo
!
! Check for those quantities for which we want yz-averages.
!
do inamex=1,nnamex
call parse_name(inamex,cnamex(inamex),cformx(inamex),'fviscmx',idiag_fviscmx)
call parse_name(inamex,cnamex(inamex),cformx(inamex),'numx',idiag_numx)
enddo
!
! Check for those quantities for which we want z-averages
!
do ixy=1,nnamexy
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'fviscmxy',idiag_fviscmxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'fviscymxy',idiag_fviscymxy)
call parse_name(ixy,cnamexy(ixy),cformxy(ixy),'fviscsmmxy',idiag_fviscsmmxy)
enddo
! check for those quantities for which we want phi-averages
!
do irz=1,nnamerz
call parse_name(irz,cnamerz(irz),cformrz(irz),'fviscrsphmphi',idiag_fviscrsphmphi)
enddo
!
call keep_compiler_quiet(lwrite)
!
endsubroutine rprint_viscosity
!***********************************************************************
subroutine pencil_criteria_viscosity
!
! All pencils that the Viscosity module depends on are specified here.
!
! 20-11-04/anders: coded
! 18-05-12/MR: request of sij2 for lvisc_simplified.and.lboussinesq added
! of graddivu for lboussinesq diasabled
!
if ((lentropy.or.ltemperature) .and. &
(lvisc_rho_nu_const .or. lvisc_rho_nu_const_bulk .or. &
lvisc_sqrtrho_nu_const .or. lvisc_nu_cspeed .or.&
lvisc_nu_const .or. lvisc_nu_tdep .or. lvisc_nu_shock .or. &
lvisc_nu_prof .or. lvisc_nu_profx .or. lvisc_spitzer .or. &
lvisc_nu_profr .or. lvisc_nu_profr_powerlaw .or. lvisc_nu_profy_bound .or. &
lvisc_nu_profr_twosteps .or. lvisc_nu_shock_profz .or. &
lvisc_nut_from_magnetic .or. lvisc_nu_shock_profr .or. lvisc_mu_cspeed))&
lpenc_requested(i_TT1)=.true.
if (lvisc_rho_nu_const .or. lvisc_rho_nu_const_bulk .or. &
lvisc_sqrtrho_nu_const .or. &
lvisc_nu_const .or. lvisc_nu_tdep .or. lvisc_nu_cspeed .or. &
lvisc_nu_prof.or.lvisc_nu_profx.or.lvisc_spitzer .or. &
lvisc_nu_profr.or.lvisc_nu_profr_powerlaw .or. &
lvisc_nu_profr_twosteps .or. lvisc_nu_profy_bound .or. &
lvisc_nut_from_magnetic.or.lvisc_mu_cspeed.or. &
(lvisc_simplified.and.lboussinesq) ) then
if ((lenergy.and.lviscosity_heat) .or. &
idiag_epsK/=0 .or. idiag_epsKint/=0 .or. idiag_epsK_LES/=0 .or. &
idiag_epsKmz/=0 .or. &
idiag_fviscmz/=0.or.idiag_fviscsmmz/=0.or.idiag_fviscmx/=0) &
lpenc_requested(i_sij2)=.true.
lpenc_requested(i_graddivu)=.true.
endif
if (lthermal_energy.and.lviscosity_heat) lpenc_requested(i_rho)=.true.
if ((lentropy.or.ltemperature).and. &
(lvisc_nu_cspeed.or.lvisc_mu_cspeed.or.lvisc_spitzer)) &
lpenc_requested(i_lnTT)=.true.
if (lvisc_smag .or.lvisc_smag_simplified .or. lvisc_smag_cross_simplified) &
lpenc_requested(i_graddivu)=.true.
if (lvisc_smag .or. lvisc_smag_simplified) lpenc_requested(i_sij2)=.true.
if (lvisc_smag_cross_simplified) lpenc_requested(i_ss12)=.true.
if (lvisc_nu_prof) lpenc_requested(i_z_mn)=.true.
if (lvisc_nu_profx) lpenc_requested(i_x_mn)=.true.
if (lvisc_nu_profr .or. lvisc_nu_profr_twosteps) then
if (lsphere_in_a_box.or.lspherical_coords) then
lpenc_requested(i_r_mn)=.true.
if (lsphere_in_a_box) lpenc_requested(i_r_mn1)=.true.
else
lpenc_requested(i_rcyl_mn)=.true.
if (lcylinder_in_a_box) lpenc_requested(i_rcyl_mn1)=.true.
endif
endif
if (lvisc_nu_profy_bound) then
if (lspherical_coords) lpenc_requested(i_r_mn)=.true.
if (lcylindrical_coords) lpenc_requested(i_rcyl_mn)=.true.
endif
if (lvisc_nu_non_newtonian) then
lpenc_requested(i_nu)=.true.
lpenc_requested(i_del2u)=.true.
lpenc_requested(i_sij)=.true.
lpenc_requested(i_sij2)=.true.
lpenc_requested(i_uijk)=.true.
endif
if (lvisc_nu_profr_powerlaw) lpenc_requested(i_rcyl_mn)=.true.
if (lvisc_nu_profr_powerlaw.and.luse_nu_rmn_prof) &
lpenc_requested(i_r_mn)=.true.
if (lvisc_rho_nu_const .or. &
lvisc_sqrtrho_nu_const .or. lvisc_nu_const .or. lvisc_nu_tdep .or. &
lvisc_smag .or. lvisc_smag_simplified .or. lvisc_smag_cross_simplified .or. &
lvisc_nu_prof .or. lvisc_nu_profx .or. lvisc_spitzer .or. &
lvisc_nu_profr_powerlaw .or. lvisc_nu_profr .or. &
lvisc_nu_profr_twosteps .or. lvisc_nu_profy_bound .or. &
lvisc_nut_from_magnetic .or. lvisc_nu_cspeed .or. lvisc_mu_cspeed) &
lpenc_requested(i_del2u)=.true.
if (.not. limplicit_viscosity .and. lvisc_simplified) lpenc_requested(i_del2u)=.true.
if (lvisc_hyper3_simplified .or. lvisc_hyper3_simplified_tdep .or. &
lvisc_hyper3_rho_nu_const .or. &
lvisc_hyper3_nu_const .or. lvisc_hyper3_rho_nu_const_symm ) &
lpenc_requested(i_del6u)=.true.
if (lvisc_hyper3_rho_nu_const_symm) then
lpenc_requested(i_grad5divu)=.true.
if ((lentropy.or.ltemperature).and.lviscosity_heat) then
lpenc_requested(i_uij5)=.true.
lpenc_requested(i_uij)=.true.
endif
endif
if (lvisc_hyper3_rho_nu_const_bulk) lpenc_requested(i_del6u_bulk)=.true.
if (lvisc_hyper2_simplified .or. lvisc_hyper2_simplified_tdep) lpenc_requested(i_del4u)=.true.
if (lvisc_rho_nu_const .or. lvisc_rho_nu_const_bulk .or. &
lvisc_sqrtrho_nu_const .or. &
lvisc_hyper3_rho_nu_const .or. lvisc_nu_cspeed .or.&
lvisc_hyper3_rho_nu_const_bulk .or. &
lvisc_hyper3_rho_nu_const_aniso .or. &
lvisc_smag .or.lvisc_smag_simplified .or. lvisc_smag_cross_simplified .or. &
lvisc_hyper3_rho_nu_const_symm .or. &
lvisc_hyper3_mu_const_strict .or. lvisc_mu_cspeed .or. &
lvisc_spitzer .or. lvisc_hyper3_cmu_const_strt_otf) &
lpenc_requested(i_rho1)=.true.
if (lvisc_hyper3_csmesh) lpenc_requested(i_cs2)=.true.
!
if (lvisc_slope_limited) then
if (lviscosity_heat) lpenc_requested(i_rho)=.true.
endif
!
if (lvisc_hyper3_cmu_const_strt_otf) then
lpenc_requested(i_del6u_strict)=.true.
lpenc_requested(i_del4graddivu)=.true.
endif
!
if (lvisc_nu_const .or. lvisc_nu_tdep .or. &
lvisc_nu_prof .or. lvisc_nu_profx .or. lvisc_nu_profy_bound .or. &
lvisc_smag .or. lvisc_smag_simplified .or. lvisc_smag_cross_simplified .or. &
lvisc_nu_profr_powerlaw .or. lvisc_nu_profr .or. &
lvisc_nu_profr_twosteps .or. lvisc_sqrtrho_nu_const .or. &
lvisc_nut_from_magnetic.or.lvisc_nu_cspeed) &
lpenc_requested(i_sglnrho)=.true.
!
if (lvisc_smag_Ma) lpenc_requested(i_Ma2)=.true.
!
if (lvisc_spitzer .or. lvisc_mu_cspeed .or. lvisc_nu_cspeed) &
lpenc_requested(i_sglnTT)=.true.
!
if (lvisc_nu_const .and. lmagfield_nu) lpenc_requested(i_b2)=.true.
if (lvisc_hyper3_nu_const) lpenc_requested(i_uij5glnrho)=.true.
if (ldensity.and.lvisc_nu_shock) then
lpenc_requested(i_graddivu)=.true.
lpenc_requested(i_shock)=.true.
lpenc_requested(i_gshock)=.true.
lpenc_requested(i_divu)=.true.
lpenc_requested(i_glnrho)=.true.
endif
if (ldensity.and.lvisc_nu_shock_profz) then
lpenc_requested(i_graddivu)=.true.
lpenc_requested(i_shock)=.true.
lpenc_requested(i_gshock)=.true.
lpenc_requested(i_divu)=.true.
lpenc_requested(i_glnrho)=.true.
endif
if (ldensity.and.lvisc_nu_shock_profr) then
lpenc_requested(i_graddivu)=.true.
lpenc_requested(i_shock)=.true.
lpenc_requested(i_gshock)=.true.
lpenc_requested(i_divu)=.true.
lpenc_requested(i_glnrho)=.true.
endif
shksmp: if (lvisc_shock_simple) then
lpenc_requested(i_shock) = .true.
lpenc_requested(i_gshock) = .true.
lpenc_requested(i_uij) = .true.
lpenc_requested(i_del2u) = .true.
endif shksmp
if (llambda_effect) then
lpenc_requested(i_uij)=.true.
lpenc_requested(i_glnrho)=.true.
!if (lKit_Olem) lpenc_requested(i_dsdr) = .true.
!if (lKit_Olem) lpenc_requested(i_gss) = .true.
endif
if (lvisc_mixture) then
lpenc_requested(i_graddivu)=.true.
lpenc_requested(i_del2u)=.true.
lpenc_requested(i_glnrho)=.true.
lpenc_requested(i_sglnrho)=.true.
lpenc_requested(i_nu)=.true.
lpenc_requested(i_gradnu)=.true.
lpenc_requested(i_sij)=.true.
if ((lentropy.or.ltemperature).and.lviscosity_heat) &
lpenc_requested(i_sij2)=.true.
endif
!
if (idiag_meshRemax/=0.or.idiag_Reshock/=0) lpenc_diagnos(i_u2)=.true.
if (idiag_visc_heatm/=0) then
lpenc_diagnos(i_visc_heat)=.true.
! lpenc_diagnos(i_sij2)=.true.
endif
if (idiag_epsK/=0 .or. idiag_epsKint/=0 .or. idiag_epsK_LES/=0 .or. &
idiag_epsKmz/=0 .or. &
idiag_viscforcezmz/=0.or.idiag_viscforcezupmz/=0.or. &
idiag_viscforcezdownmz/=0) then
lpenc_diagnos(i_rho)=.true.
! lpenc_diagnos(i_sij2)=.true.
! JW: There are also viscosities, which do not need sij
! PJK: But what if you don't use one of those?
! JW: then they are anyway requested.
endif
if (idiag_sijoiojm/=0) then
lpenc_diagnos(i_oo)=.true.
lpenc_diagnos(i_sij)=.true.
endif
if (idiag_Sij2m/=0.) lpenc_diagnos(i_sij2)=.true.
if (idiag_epsK/=0 .or. idiag_epsKint/=0 .or. idiag_epsKmz/=0 .or. idiag_epsK_LES/=0) then
lpenc_diagnos(i_visc_heat)=.true.
lpenc_diagnos(i_uu)=.true.
endif
if (lvisc_nu_shock.and.(idiag_epsK/=0 .or. idiag_epsKint/=0)) then
lpenc_diagnos(i_fvisc)=.true.
lpenc_diagnos(i_diffus_total)=.true.
lpenc_diagnos(i_shock)=.true.
lpenc_diagnos(i_divu)=.true.
lpenc_diagnos(i_rho)=.true.
endif
if (lvisc_nu_shock_profz.and.(idiag_epsK/=0 .or. idiag_epsKint/=0)) then
lpenc_diagnos(i_fvisc)=.true.
lpenc_diagnos(i_diffus_total)=.true.
lpenc_diagnos(i_shock)=.true.
lpenc_diagnos(i_divu)=.true.
lpenc_diagnos(i_rho)=.true.
endif
if (lvisc_nu_shock_profr.and.(idiag_epsK/=0 .or. idiag_epsKint/=0)) then
lpenc_diagnos(i_fvisc)=.true.
lpenc_diagnos(i_diffus_total)=.true.
lpenc_diagnos(i_shock)=.true.
lpenc_diagnos(i_divu)=.true.
lpenc_diagnos(i_rho)=.true.
endif
if (lvisc_heat_as_aux) then
lpenc_diagnos(i_visc_heat)=.true.
! lpenc_diagnos(i_rho)=.true.
! lpenc_diagnos(i_sij2)=.true.
endif
if ((idiag_meshRemax/=0.or.idiag_dtnu/=0).and.(lvisc_nu_shock &
.or.lvisc_nu_shock_profz.or.lvisc_nu_shock_profr)) then
lpenc_diagnos(i_shock)=.true.
endif
if (idiag_qfviscm/=0) then
lpenc_diagnos(i_curlo)=.true.
lpenc_diagnos(i_fvisc)=.true.
endif
if (idiag_ufviscm/=0) then
lpenc_diagnos(i_uu)=.true.
lpenc_diagnos(i_fvisc)=.true.
endif
if (idiag_Reshock/=0) lpenc_diagnos(i_shock)=.true.
if (idiag_fviscmz/=0.or.idiag_fviscsmmz/=0.or.idiag_fviscmx/=0) then
lpenc_diagnos(i_uu)=.true.
! lpenc_diagnos(i_rho)=.true.
! lpenc_diagnos(i_sij)=.true.
endif
! JW: There are also viscosities, which do not need sij or rho
! PJK: But what if you don't use one of those?
! JW: p%uu is needed, these others are already requested.
if (idiag_numx/=0) then
lpenc_diagnos(i_nu) = .true.
endif
if (idiag_fviscmxy/=0 .or. idiag_fviscymxy/=0 .or. &
idiag_fviscsmmxy/=0) then
! lpenc_diagnos2d(i_nu_smag)=.true.
lpenc_diagnos2d(i_uu)=.true.
! lpenc_diagnos2d(i_rho)=.true.
! lpenc_diagnos2d(i_sij)=.true.
endif
! JW: There are also viscosities, which do not need sij or rho
! PJK: But what if you don't use one of those?
! JW: p%uu is needed, these others are already requested.
if (idiag_fviscrsphmphi/=0) then
lpenc_diagnos2d(i_evr)=.true.
endif
if (lboussinesq) lpenc_requested(i_graddivu)=.false.
if (damp_sound/=0.) lpenc_requested(i_divu)=.true.
if (lvisc_hyper3_mesh_residual) lpenc_requested(i_der6u_res)=.true.
!
endsubroutine pencil_criteria_viscosity
!***********************************************************************
subroutine pencil_interdep_viscosity(lpencil_in)
!
! Interdependency among pencils from the Viscosity module is specified here.
!
! 20-11-04/anders: coded
!
logical, dimension (npencils) :: lpencil_in
!
if (lvisc_simplified .and. lpencil_in(i_visc_heat)) lpencil_in(i_o2)=.true.
!
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.
!
! 20-nov-04/anders: coded
! 18-may-12/MR: calculation of viscous heat for boussinesq added
! 14-oct-15/MR: corrected viscous force for slope-limited flux
! 03-apr-20/joern: restructured and fixed slope-limited diffusion
!
use Deriv, only: der5i1j,der6
use Diagnostics, only: max_mn_name, sum_mn_name
use Sub
!
real, dimension (mx,my,mz,mfarray) :: f
type (pencil_case) :: p
intent(inout) :: f,p
!
real, dimension (nx,3) :: tmp,tmp2,gradnu,sgradnu,gradnu_shock
real, dimension (nx) :: murho1,zetarho1,muTT,tmp3,tmp4,pnu_shock
real, dimension (nx) :: lambda_phi,prof,prof2,derprof,derprof2,qfvisc
real, dimension (nx) :: gradnu_effective,fac,advec_hypermesh_uu
real, dimension (nx,3) :: deljskl2,fvisc_nnewton2
real, dimension (nx,3,3) :: d_sld_flux
!
integer :: i,j,ju,ii,jj,kk,ll
!
! Viscous force and viscous heating are calculated here (for later use).
!
p%fvisc=0.0 !!! not needed
if (lpencil(i_visc_heat)) p%visc_heat=0.0
if (lfirst .and. ldt) then
p%diffus_total=0.0
p%diffus_total2=0.0
p%diffus_total3=0.0
endif
!
! viscous force: nu*del2v
! -- not physically correct (no momentum conservation), but
! numerically easy and in most cases qualitatively OK.
! For boussinesq (divu=0) this is exact.
! In all other cases we use nu*o2.
!
if (.not. limplicit_viscosity .and. lvisc_simplified) then
p%fvisc=p%fvisc+nu*p%del2u
if (lpencil(i_visc_heat)) then
if (lboussinesq) then
p%visc_heat=p%visc_heat+2.*nu*p%sij2
else
p%visc_heat=p%visc_heat+nu*p%o2
endif
endif
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+nu
endif
!
! Non-newtonian viscosity
!
if (lvisc_nu_non_newtonian) then
! print*,'viscous force: div(nu*Sij); with'
! print*,'nu a function of square of strain rate'
! fvisc_i = \partial_j \nu(S_{kl}S_{kl}) S_{ij}
! = \nu(S_{kl}S_{kl}) \partial_j S_{ij} + S_{ij}\partial_j \nu(S_{kl}S_{kl})
! = \nu(S_{kl}S_{kl}) \partial_j (uij+uji) + S_{ij}\nu^{\prime} \partial_j [S_{kl}S_{kl}]
! = \nu(S_{kl}S_{kl}) \partial_jj (ui) + S_{ij}\nu^{\prime} \partial_j [(ukl+ulk)(ukl+ulk)]
! = \nu(.) \del2 (ui) + S_{ij}\nu^{\prime} \partial_j [ukl ukl +ulk ukl + ukl ulk + ulk ulk]
! = \nu(.) \del2 (ui) + S_{ij}\nu^{\prime} [(uklj ukl + ukl uklj) + (ulkj ukl + ulk uklj)
! +(uklj ulk + ukl ulkj) + (ulkj ulk + ulk ulkj)]
! = \nu(.) \del2 (ui) + S_{ij}\nu^{\prime} 2[uklj S_{lk}+ ulkj S_{kl}]
!
deljskl2=0.
gradnu_effective=0.
do jj=1,3
do kk=1,3; do ll=1,3
deljskl2 (:,jj) = deljskl2 (:,jj) + &
p%uijk(:,kk,ll,jj)*p%sij(:,ll,kk) + p%uijk(:,ll,kk,jj)*p%sij(:,kk,ll)
enddo;enddo
enddo
do ii=1,3
fvisc_nnewton2(:,ii) = sum(p%sij(:,ii,:)*deljskl2,2)
enddo
call getnu_non_newtonian(p%sij2,p%nu,gradnu_effective)
do ii=1,3
p%fvisc(:,ii)=p%nu(:)*p%del2u(:,ii) + gradnu_effective(:)*fvisc_nnewton2(:,ii)
enddo
!
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+nu
endif
!
! viscous force: mu/rho*(del2u+graddivu/3)
! -- the correct expression for rho*nu=const
! As a test for vorticity generation, we also allow for the possibility
! of setting the prefactor to mu (without 1/rho factor). In that case
! we set lvisc_rho_nu_const_prefact=T
!
if (lvisc_rho_nu_const) then
if (lvisc_rho_nu_const_prefact) then
murho1=nu !(=mu=dynamical viscosity)
else
murho1=nu*p%rho1 !(=mu/rho)
endif
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i) + &
murho1*(p%del2u(:,i)+1.0/3.0*p%graddivu(:,i))
enddo
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*murho1*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+murho1
endif
!
! viscous force: zeta/rho*graddivu (constant dynamic bulk viscosity)
!
if (lvisc_rho_nu_const_bulk) then
zetarho1=zeta*p%rho1 !(=zeta/rho)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+zetarho1*p%graddivu(:,i)
enddo
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+zetarho1*p%divu**2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+zetarho1
endif
!
! viscous force: mu/sqrt(rho)*(del2u+graddivu/3)
! [Implemented by Fred Gent in r13626 of 13 April 2010]
! [AB: this may not correspond to a symmetric stress tensor]
! PJK: 21-09-13 modified to include the missing term
! viscous force: mu/sqrt(rho)*(del2u+graddivu/3+S.glnrho)
!
if (lvisc_sqrtrho_nu_const) then
murho1=nu*sqrt(p%rho1)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i) + &
murho1*(p%del2u(:,i)+1.0/3.0*p%graddivu(:,i) &
+ p%sglnrho(:,i))
enddo
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*murho1*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+murho1
endif
!
! viscous force: nu*sqrt(TT)/rho*(del2u+graddivu/3+S.glnTT)
! fred: 23.9.17 replaced 0.5 with nu_cspeed so exponent can be generalised
!
if (lvisc_mu_cspeed) then
muTT=nu*p%rho1*exp(nu_cspeed*p%lnTT)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i) + &
muTT*(p%del2u(:,i)+1.0/3.0*p%graddivu(:,i)&
+2*nu_cspeed*p%sglnTT(:,i))
enddo
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*muTT*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+muTT
endif
!
! viscous force: nu*TT^2.5/rho*(del2u+graddivu/3+5S.glnTT)
!
if (lvisc_spitzer) then
if (nu_spitzer_max .ne. 0.0) then
tmp3=nu_spitzer*p%rho1*exp(2.5*p%lnTT)
muTT=nu_spitzer_max*tmp3/sqrt(nu_spitzer_max**2.+tmp3**2.)
else
muTT=nu_spitzer*p%rho1*exp(2.5*p%lnTT)
endif
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i) + &
muTT*(p%del2u(:,i)+1.0/3.0*p%graddivu(:,i)+5.*p%sglnTT(:,i))
enddo
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*muTT*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+muTT
endif
!
! viscous force: nu*sqrt(TT)*(del2u+graddivu/3+2S.glnrho)
! -- for numerical stability viscous force propto soundspeed in interstellar run
! fred: 23.9.17 replaced 0.5 with nu_cspeed so exponent can be generalised
!
if (lvisc_nu_cspeed) then
muTT=nu*exp(nu_cspeed*p%lnTT)
if (ldensity) then
do i=1,3
p%fvisc(:,i) = p%fvisc(:,i) + 2*muTT*p%sglnrho(:,i)&
+muTT*(p%del2u(:,i) + 1./3.*p%graddivu(:,i)&
+2*nu_cspeed*p%sglnTT(:,i))
enddo
! Tobi: This is not quite the full story in the presence of linear
! shear. In this case the rate-of-strain tensor S has xy and yx
! components
! S_xy = S_yx = 1/2 (du_x/dy + du_y/dx + Sshear)
! so that one must add
!if (.false.) then
! p%fvisc(:,1) = p%fvisc(:,1) + Sshear*muTT*p%glnrho(:,2)
! p%fvisc(:,2) = p%fvisc(:,2) + Sshear*muTT*p%glnrho(:,1)
!endif
else
if (lmeanfield_nu) then
if (meanfield_nuB/=0.) then
call multsv_mn(muTT,p%del2u+1./3.*p%graddivu,tmp)
call multsv_mn_add(1./sqrt(1.+p%b2/meanfield_nuB**2),&
p%fvisc+tmp,p%fvisc)
endif
else
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+muTT*(p%del2u(:,i)&
+1.0/3.0*p%graddivu(:,i)+p%sglnTT(:,i))
enddo
endif
endif
!
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*muTT*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+muTT
endif
!
! viscous force: nu*(del2u+graddivu/3+2S.glnrho)
! -- the correct expression for nu=const
!
if (lvisc_nu_const) then
if (ldensity) then
fac=nu
if (damp_sound/=0.) fac = fac+damp_sound*abs(p%divu)
do j=1,3
p%fvisc(:,j) = p%fvisc(:,j) + fac*(p%del2u(:,j) + 2.*p%sglnrho(:,j) + 1./3.*p%graddivu(:,j))
enddo
! Tobi: This is not quite the full story in the presence of linear
! shear. In this case the rate-of-strain tensor S has xy and yx
! components
! S_xy = S_yx = 1/2 (du_x/dy + du_y/dx + Sshear)
! so that one must add
!if (.false.) then
! p%fvisc(:,1) = p%fvisc(:,1) + Sshear*nu*p%glnrho(:,2)
! p%fvisc(:,2) = p%fvisc(:,2) + Sshear*nu*p%glnrho(:,1)
!endif
else
if (lmeanfield_nu) then
if (meanfield_nuB/=0.) then
call multsv_mn_add(1./sqrt(1.+p%b2/meanfield_nuB**2), &
p%fvisc+nu*(p%del2u+1./3.*p%graddivu),p%fvisc)
endif
else
p%fvisc=p%fvisc+nu*(p%del2u+1.0/3.0*p%graddivu)
endif
endif
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*nu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+nu
endif
!
! Viscous force: nu(t)*(del2u+graddivu/3+2S.glnrho) [correct for nu=const].
! The following allows us to let nu change with time, t-nu_tdep_toffset.
! The nu_tdep_toffset is used in cosmology where time starts at t=1.
! lresi_nu_tdep_t0_norm is not the default because of backward compatbility.
! The default is problematic because then nu_tdep /= nu for t < nu_tdep_t0.
!
if (lvisc_nu_tdep .or. lvisc_hyper3_simplified_tdep) then
if (lvisc_nu_tdep_t0_norm) then
nu_tdep=nu*max(real(t-nu_tdep_toffset)/nu_tdep_t0,1.)**nu_tdep_exponent
else
nu_tdep=nu*max(real(t-nu_tdep_toffset),nu_tdep_t0)**nu_tdep_exponent
endif
p%fvisc=p%fvisc+2*nu_tdep*p%sglnrho+nu_tdep*(p%del2u+1./3.*p%graddivu)
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*nu_tdep*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+nu_tdep
endif
!
! Viscous force: nu*(del2u+graddivu/3+2S.glnrho)+2S.gradnu
!
if (lvisc_mixture) then
call multmv(p%sij,p%gradnu,sgradnu)
do i=1,3
p%fvisc(:,i)=p%nu*(p%del2u(:,i)+1./3.*p%graddivu(:,i)) + 2*sgradnu(:,i)
enddo
if (ldensity) then
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i) + 2*p%nu*p%sglnrho(:,i)
enddo
endif
!
! Viscous heating and time step.
!
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*p%nu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+p%nu
endif
!
! Viscous force: nu(x)*(del2u+graddivu/3+2S.glnrho)+2S.gnu
! -- here the nu viscosity depends on x; nu_jump=nu2/nu1
! pnu = nu + (nu*(nu_jump-1.))*step(abs(p%x_mn),xnu,widthnu)
!
if (lvisc_nu_profx.or.lvisc_nu_profr) then
!if (lvisc_nu_profx) tmp3=p%x_mn
!AB: Petri, Matthias, or somebody using spherical coordinates;
!AB: the line above seems wrong and should simply be like so:
if (lvisc_nu_profx) tmp3=x(l1:l2)
if (lvisc_nu_profr) then
if (lspherical_coords.or.lsphere_in_a_box) then
tmp3=p%r_mn
else
tmp3=p%rcyl_mn
endif
endif
if (lvisc_nu_profx.and.lvisc_nu_profr) then
print*,'You are using both radial and horizontal '
print*,'profiles for a viscosity jump. Are you sure '
print*,'this is reasonable? Better stop and check.'
call fatal_error("calc_pencils_viscosity","")
endif
pnu = nu + (nu*(nu_jump-1.))*(step(tmp3,xnu ,widthnu) - step(tmp3,xnu2,widthnu))
tmp4 = (nu*(nu_jump-1.))*(der_step(tmp3,xnu ,widthnu)-der_step(tmp3,xnu2,widthnu))
!
! Initialize gradnu before calculating it, otherwise gfortran complains.
!
gradnu=0.
call get_gradnu(tmp4,lvisc_nu_profx,&
lvisc_nu_profr,p,gradnu)
! A routine for write_xprof should be written here
! call write_xprof('visc',pnu)
!gradnu(:,1) = (nu*(nu_jump-1.))*der_step(tmp3,xnu,widthnu)
!gradnu(:,2) = 0.
!gradnu(:,3) = 0.
call multmv(p%sij,gradnu,sgradnu)
call multsv(pnu,2*p%sglnrho+p%del2u+1./3.*p%graddivu,tmp)
!tobi: The following only works with operator overloading for pencils
! (see sub.f90). Commented out because it seems to be slower.
!p%fvisc=p%fvisc+2*pnu*p%sglnrho+pnu*(p%del2u+1./3.*p%graddivu) &
! +2*sgradnu
p%fvisc=p%fvisc+tmp+2*sgradnu
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*pnu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+pnu
endif
!
! Radial viscosity profile from power law.
! Viscous force: nu(x)*(del2u+graddivu/3+2S.glnrho)+2S.gnu
! -- here the nu viscosity depends on r; nu=nu_0*(r/r0)^(-pnlaw)
!
if (lvisc_nu_profr_powerlaw) then
if (.not.luse_nu_rmn_prof) then
pnu = nu*(p%rcyl_mn/xnu)**(-pnlaw)
!
! viscosity gradient
!
if (lcylindrical_coords) then
gradnu(:,1) = -pnlaw*nu*(p%rcyl_mn/xnu)**(-pnlaw-1)*1/xnu
gradnu(:,2) = 0.
gradnu(:,3) = 0.
elseif (lspherical_coords) then
gradnu(:,1) = -pnlaw*nu*p%rcyl_mn**(-pnlaw-1)*sinth(m)
gradnu(:,2) = -pnlaw*nu*p%rcyl_mn**(-pnlaw-1)*costh(m)
gradnu(:,3) = 0.
else
print*,'power-law viscosity only implemented '
print*,'for spherical and cylindrical coordinates'
call fatal_error("calc_pencils_viscosity","")
endif
else
pnu = nu*(p%r_mn/xnu)**(-pnlaw)
!
! viscosity gradient
!
if (lspherical_coords) then
gradnu(:,1) = -pnlaw*nu*(p%r_mn/xnu)**(-pnlaw-1)*1/xnu
gradnu(:,2) = 0.
gradnu(:,3) = 0.
else
print*,'power-law viscosity with luse_nu_rmn_prof=T only '
print*,'implemented for spherical coordinates'
call fatal_error("calc_pencils_viscosity","")
endif
endif
!
call multmv(p%sij,gradnu,sgradnu)
call multsv(pnu,2*p%sglnrho+p%del2u+1./3.*p%graddivu,tmp)
p%fvisc=p%fvisc+tmp+2*sgradnu
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*pnu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+pnu
endif
!
! Radial viscosity profile with two steps; very similar to nu_profr
! Viscous force: nu(x)*(del2u+graddivu/3+2S.glnrho)+2S.gnu
! -- here the nu viscosity depends on r; nu_jump=nu2/nu1
!
if (lvisc_nu_profr_twosteps) then
if (lspherical_coords.or.lsphere_in_a_box) then
tmp3=p%r_mn
else
tmp3=p%rcyl_mn
endif
if (lvisc_nu_profx.and.lvisc_nu_profr_twosteps) then
print*,'You are using both radial and horizontal '
print*,'profiles for a viscosity jump. Are you sure '
print*,'this is reasonable? Better stop and check.'
call fatal_error("","")
endif
prof2 = step(tmp3,xnu2,widthnu2)
prof = step(tmp3,xnu,widthnu)-prof2
derprof2 = der_step(tmp3,xnu2,widthnu2)
derprof = der_step(tmp3,xnu,widthnu)-derprof2
!
pnu = nu + (nu*(nu_jump-1.))*prof+(nu*(nu_jump2-1.))*prof2
tmp4 = nu + (nu*(nu_jump-1.))*derprof+(nu*(nu_jump2-1.))*derprof2
!
! Initialize gradnu before calculating it, otherwise gfortran complains.
!
gradnu=0.
call get_gradnu(tmp4,lvisc_nu_profx,lvisc_nu_profr_twosteps,p,gradnu)
call multmv(p%sij,gradnu,sgradnu)
call multsv(pnu,2*p%sglnrho+p%del2u+1./3.*p%graddivu,tmp)
p%fvisc=p%fvisc+tmp+2*sgradnu
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*pnu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+pnu
endif
!
! horizontal viscosity profile with two steps
! enhancement at the y boundaries; very similar to nu_profr_twosteps
! dynu define how wide the enhanced viscosity layer is.
! Viscous force: nu(y)*(del2u+graddivu/3+2S.glnrho)+2S.gnu
! -- here the nu viscosity depends on y; nu_jump=nu2/nu1
!
if (lvisc_nu_profy_bound) then
if (lspherical_coords.or.lsphere_in_a_box) then
tmp3=p%r_mn
else if (lcylindrical_coords) then
tmp3=p%rcyl_mn
else
tmp3=1.0
endif
!
tmp4 = spread(y(m),1,nx)
prof = step(tmp4,xyz1(2)-3*dynu, dynu) + step(tmp4,xyz0(2)+3*dynu, -dynu)
derprof = der_step(tmp4,xyz1(2)-3*dynu, dynu) + der_step(tmp4,xyz0(2)+3*dynu, -dynu)
!
pnu = nu + (nu*(nu_jump-1.))*prof
!
gradnu(:,1) = 0.
gradnu(:,2) = (nu*(nu_jump-1.))*derprof/tmp3
gradnu(:,3) = 0.
call multmv(p%sij,gradnu,sgradnu)
call multsv(pnu,2*p%sglnrho+p%del2u+1./3.*p%graddivu,tmp)
p%fvisc=p%fvisc+tmp+2*sgradnu
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*pnu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+pnu
endif
!
! turbulent viscosity profile from magnetic
! viscous force: nu(z)*(del2u+graddivu/3+2S.glnrho)+2S.gnu
! The following can in future be generalized to read
!
if (lvisc_nut_from_magnetic) then
pnu=PrM_turb*etat_x*etat_y(m)*etat_z(n)
gradnu(:,1)=PrM_turb*detat_x*etat_y(m)*etat_z(n)
gradnu(:,2)=PrM_turb*etat_x*detat_y(m)*etat_z(n)
gradnu(:,3)=PrM_turb*etat_x*etat_y(m)*detat_z(n)
call multmv(p%sij,gradnu,sgradnu)
call multsv(pnu,2*p%sglnrho+p%del2u+1./3.*p%graddivu,tmp)
p%fvisc=p%fvisc+tmp+2*sgradnu
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*pnu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+pnu
endif
!
! viscous force: nu(z)*(del2u+graddivu/3+2S.glnrho)+2S.gnu
! -- here the nu viscosity depends on z; nu_jump=nu2/nu1
! When widthnu is negative, the viscosity increases with height.
! Otherwise, for positive widthnu, it *decreases* with height.
!
if (lvisc_nu_prof) then
pnu = nu + (nu*(nu_jump-1.))*step(p%z_mn,znu,-widthnu)
!
! This indicates that the profile is meant to be applied in Cartesian
! coordinates only. Why then z taken from pencil?
!
gradnu(:,1) = 0.
gradnu(:,2) = 0.
gradnu(:,3) = (nu*(nu_jump-1.))*der_step(p%z_mn,znu,-widthnu)
call multmv(p%sij,gradnu,sgradnu)
call multsv(pnu,2*p%sglnrho+p%del2u+1./3.*p%graddivu,tmp)
!tobi: The following only works with operator overloading for pencils
! (see sub.f90). Commented out because it seems to be slower.
!p%fvisc=p%fvisc+2*pnu*p%sglnrho+pnu*(p%del2u+1./3.*p%graddivu) &
! +2*sgradnu
p%fvisc=p%fvisc+tmp+2*sgradnu
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*pnu*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+pnu
endif
!
! viscous force: nu_shock
!
if (lvisc_nu_shock) then
if (ldensity) then
!tobi: The following only works with operator overloading for pencils
! (see sub.f90). Commented out because it seems to be slower.
!tmp=nu_shock*(p%shock*(p%divu*p%glnrho+p%graddivu)+p%divu*p%gshock)
call multsv(p%divu,p%glnrho,tmp2)
tmp=tmp2 + p%graddivu
call multsv(nu_shock*p%shock,tmp,tmp2)
call multsv_add(tmp2,nu_shock*p%divu,p%gshock,tmp)
p%fvisc=p%fvisc+tmp
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+(nu_shock*p%shock)
if (lpencil(i_visc_heat)) &
p%visc_heat=p%visc_heat+nu_shock*p%shock*p%divu**2
endif
endif
!
! viscous force: nu_shock with vertical profile
!
if (lvisc_nu_shock_profz) then
if (ldensity) then
pnu_shock = nu_shock + (nu_shock*(nu_jump_shock-1.))*step(p%z_mn,znu_shock,-widthnu_shock)
!
! This indicates that the profile is meant to be applied in Cartesian
! coordinates only. Why then z taken from pencil?
!
gradnu_shock(:,1) = 0.
gradnu_shock(:,2) = 0.
gradnu_shock(:,3) = (nu_shock*(nu_jump_shock-1.))*der_step(p%z_mn,znu_shock,-widthnu_shock)
!
call multsv(p%divu,p%glnrho,tmp2)
tmp=tmp2 + p%graddivu
call multsv(pnu_shock*p%shock,tmp,tmp2)
call multsv_add(tmp2,pnu_shock*p%divu,p%gshock,tmp)
call multsv_mn_add(p%shock*p%divu,gradnu_shock,tmp)
p%fvisc=p%fvisc+tmp
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+(pnu_shock*p%shock)
if (lpencil(i_visc_heat)) &
p%visc_heat=p%visc_heat+pnu_shock*p%shock*p%divu**2
endif
endif
!
! viscous force: nu_shock with radial profile
!
if (lvisc_nu_shock_profr) then
if (ldensity) then
if (lspherical_coords.or.lsphere_in_a_box) then
tmp3=p%r_mn
else
tmp3=p%rcyl_mn
endif
pnu_shock = nu_shock + (nu_shock*(nu_jump_shock-1.))*step(tmp3,xnu_shock,widthnu_shock)
!
! This indicates that the profile is meant to be applied in spherical and
! cylindrical coordinates only, referring to r and pomega, respectively.
! Hence not correct for 'sphere in a box'!
! Why then r taken from pencil?
!
gradnu_shock(:,1) = (nu_shock*(nu_jump_shock-1.))*der_step(tmp3,xnu_shock,widthnu_shock)
gradnu_shock(:,2) = 0.
gradnu_shock(:,3) = 0.
!
call multsv(p%divu,p%glnrho,tmp2)
tmp=tmp2 + p%graddivu
call multsv(pnu_shock*p%shock,tmp,tmp2)
call multsv_add(tmp2,pnu_shock*p%divu,p%gshock,tmp)
call multsv_mn_add(p%shock*p%divu,gradnu_shock,tmp)
p%fvisc=p%fvisc+tmp
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+(pnu_shock*p%shock)
if (lpencil(i_visc_heat)) &
p%visc_heat=p%visc_heat+pnu_shock*p%shock*p%divu**2
endif
endif
!
! viscous force: div(nu_shock * grad(uu_i))
!
shksmp: if (lvisc_shock_simple) then
do i = 1, 3
call dot(p%gshock, p%uij(:,i,:), tmp3)
tmp(:,i) = tmp3 + p%shock * p%del2u(:,i)
enddo
p%fvisc = p%fvisc + nu_shock * tmp
if (lfirst .and. ldt) p%diffus_total = p%diffus_total + nu_shock * p%shock
if (lpencil(i_visc_heat) .and. headtt) call warning('calc_pencils_viscosity', 'shock heating does not exist. ')
endif shksmp
!
! viscous force: nu_hyper2*de46v (not momentum-conserving)
!
if (lvisc_hyper2_simplified) then
p%fvisc=p%fvisc-nu_hyper2*p%del4u
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper2_simplified')
endif
endif
if (lfirst .and. ldt) p%diffus_total2=p%diffus_total2+nu_hyper2
endif
!
! viscous force: nu_hyper3*del6v (not momentum-conserving)
!
if (lvisc_hyper3_simplified) then
p%fvisc=p%fvisc+nu_hyper3*p%del6u
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_simplified')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_hyper3
endif
!
! viscous force with t-dependent nu: nu*del6v (not momentum-conserving)
!
if (lvisc_hyper2_simplified_tdep) then
p%fvisc=p%fvisc-nu_tdep*p%del4u
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper2_simplified')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_tdep
endif
!
if (lvisc_hyper3_simplified_tdep) then
p%fvisc=p%fvisc+nu_tdep*p%del6u
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_simplified')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_tdep
endif
!
! General way of coding an anisotropic hyperviscosity.
!
if (lvisc_hyper3_polar) then
do j=1,3
ju=j+iuu-1
do i=1,3
call der6(f,ju,tmp3,i,IGNOREDX=.true.)
p%fvisc(:,j) = p%fvisc(:,j) + &
nu_hyper3*pi4_1*tmp3*dline_1(:,i)**2
enddo
enddo
if (lpencil(i_visc_heat)) then
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '//&
'is not implemented for lvisc_hyper3_polar')
endif
endif
if (lfirst .and. ldt) &
p%diffus_total3=p%diffus_total3+nu_hyper3*pi4_1*dxmin_pencil**4
endif
!
! Following Axel's hyper3_mesh for density
!
if (lvisc_hyper3_mesh) then
do j=1,3
ju=j+iuu-1
do i=1,3
call der6(f,ju,tmp3,i,IGNOREDX=.true.)
if (ldynamical_diffusion) then
p%fvisc(:,j) = p%fvisc(:,j) + nu_hyper3_mesh * tmp3 * dline_1(:,i)
else
p%fvisc(:,j) = p%fvisc(:,j) + &
nu_hyper3_mesh*pi5_1/60.* tmp3 *dline_1(:,i)
endif
enddo
enddo
if (lpencil(i_visc_heat)) then
if (headtt) &
call warning('calc_pencils_viscosity', 'viscous heating term '//&
'is not implemented for lvisc_hyper3_mesh')
endif
if (lfirst .and. ldt) then
if (ldynamical_diffusion) then
p%diffus_total3 = p%diffus_total3 + nu_hyper3_mesh
advec_hypermesh_uu = 0.0
else
advec_hypermesh_uu=nu_hyper3_mesh*pi5_1*sqrt(dxyz_2)
endif
advec2_hypermesh=advec2_hypermesh+advec_hypermesh_uu**2
endif
endif
!
! Hyper-mesh diffusion with residual velocity instead of the full velocity.
!
if (lvisc_hyper3_mesh_residual) then
do j=1,3
ju=j+iuu-1
if (ldynamical_diffusion) then
p%fvisc(:,j) = p%fvisc(:,j) + nu_hyper3_mesh * sum(p%der6u_res(:,:,j) * dline_1,2)
else
p%fvisc(:,j) = p%fvisc(:,j) + &
nu_hyper3_mesh*pi5_1/60.*sum(p%der6u_res(:,:,j)*dline_1,2)
endif
enddo
if (lpencil(i_visc_heat)) then
if (headtt) &
call warning('calc_pencils_viscosity', 'viscous heating term '//&
'is not implemented for lvisc_hyper3_mesh')
endif
if (lfirst .and. ldt) then
if (ldynamical_diffusion) then
p%diffus_total3 = p%diffus_total3 + nu_hyper3_mesh
advec_hypermesh_uu = 0.0
else
advec_hypermesh_uu=nu_hyper3_mesh*pi5_1*sqrt(dxyz_2)
endif
advec2_hypermesh=advec2_hypermesh+advec_hypermesh_uu**2
endif
endif
!
if (lvisc_hyper3_csmesh) then
do j=1,3
ju=j+iuu-1
do i=1,3
call der6(f,ju,tmp3,i,IGNOREDX=.true.)
if (ldynamical_diffusion) then
p%fvisc(:,j)=p%fvisc(:,j)+nu_hyper3_mesh*sqrt(p%cs2) &
*tmp3*dline_1(:,i)
else
p%fvisc(:,j)=p%fvisc(:,j)+nu_hyper3_mesh*sqrt(p%cs2) &
*pi5_1/60.*tmp3*dline_1(:,i)
endif
enddo
enddo
if (lpencil(i_visc_heat)) then
if (headtt) &
call warning('calc_pencils_viscosity', 'viscous heating term '//&
'is not implemented for lvisc_hyper3_csmesh')
endif
if (lfirst .and. ldt) then
if (ldynamical_diffusion) then
p%diffus_total3=p%diffus_total3+nu_hyper3_mesh*sqrt(p%cs2)
advec_hypermesh_uu=0.0
else
advec_hypermesh_uu=nu_hyper3_mesh*pi5_1*sqrt(dxyz_2*p%cs2)
endif
advec2_hypermesh=advec2_hypermesh+advec_hypermesh_uu**2
endif
endif
!
! viscous force: mu/rho*del6u
!
if (lvisc_hyper3_rho_nu_const) then
murho1=nu_hyper3*p%rho1 ! (=mu_hyper3/rho)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+murho1*p%del6u(:,i)
enddo
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_rho_nu_const')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+murho1
endif
!
! For tau_ij=d^5u_i/dx_j^5 + d^5u_j/dx_i^5
! Viscous force: du/dt = mu/rho*{del6(u) + grad5[div(u)]}
!
if (lvisc_hyper3_rho_nu_const_symm) then
murho1=nu_hyper3*p%rho1 ! (=mu_hyper3/rho)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+murho1*(p%del6u(:,i) + p%grad5divu(:,i))
enddo
if (lpencil(i_visc_heat)) then
do i=1,3; do j=1,3
!
! Dissipation is *not* positive definite.
!
p%visc_heat=p%visc_heat + &
0.5*nu_hyper3*(p%uij5(:,i,j)+p%uij5(:,j,i))*p%uij(:,i,j)
enddo; enddo
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+murho1
endif
!
! Viscous force:
! du/dt = mu/rho*{del2(del2(del2(u))) + 1/3*del2(del2(grad(div(u))
! This viscosity type requires HYPERVISC_STRICT = hypervisc_strict_fft in
! your Makefile.local.
!
if (lvisc_hyper3_mu_const_strict) then
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+nu_hyper3*p%rho1*f(l1:l2,m,n,ihypvis-1+i)
enddo
if (lpencil(i_visc_heat)) then ! Should be eps=2*mu*{del2[del2(S)]}^2
if (headtt) then ! (see Haugen & Brandenburg 2004 eq. 7)
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_mu_const_strict')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_hyper3
endif
!
! A strict hyperviscosity, with the mixed-derivatives coded in deriv, and assessed
! on the fly. Much faster than the one above using the hypervisc_strict module
!
if (lvisc_hyper3_cmu_const_strt_otf) then
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+nu_hyper3*p%rho1*(p%del6u_strict(:,i) + 1./3*p%del4graddivu(:,i))
enddo
if (lpencil(i_visc_heat)) then ! Should be eps=2*mu*{del2[del2(S)]}^2
if (headtt) then ! (see Haugen & Brandenburg 2004 eq. 7)
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_mu_const_strict_otf')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_hyper3
endif
!
! Viscous force:
! du/dt = 1/rho*div[2*rho*nu_3*S^(3)]
! This viscosity type requires HYPERVISC_STRICT = hypervisc_strict_fft in
! your Makefile.local.
!
if (lvisc_hyper3_nu_const_strict) then
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+nu_hyper3*f(l1:l2,m,n,ihypvis-1+i)
enddo
if (lpencil(i_visc_heat)) then
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_nu_const_strict')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_hyper3
endif
!
! viscous force: f_i = mu_i/rho*del6u
! Used for non-cubic cells
!
if (lvisc_hyper3_rho_nu_const_aniso) then
call del6fjv(f,nu_aniso_hyper3,iuu,tmp)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+tmp(:,i)*p%rho1
enddo
!
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_rho_nu_const_aniso')
endif
endif
!
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+&
(nu_aniso_hyper3(1)*dline_1(:,1)**6 + &
nu_aniso_hyper3(2)*dline_1(:,2)**6 + &
nu_aniso_hyper3(3)*dline_1(:,3)**6)/ dxyz_6
!
endif
!
! viscous force: f_i = (nu_j.del6)u_i + nu_j.uij5.glnrho
! Used for non-cubic cells
!
if (lvisc_hyper3_nu_const_aniso) then
call del6fjv(f,nu_aniso_hyper3,iuu,tmp)
!
do i=1,3
tmp3=0.
do j=1,3
tmp3=tmp3+p%uij(:,i,j)*p%glnrho(:,j)*nu_aniso_hyper3(j)
enddo
!
p%fvisc(:,i)=p%fvisc(:,i)+tmp(:,i)+tmp3
enddo
!
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_nu_const_aniso')
endif
endif
!
! diffusion time: it will be multiplied by dxyz_2 again further down
!
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+&
(nu_aniso_hyper3(1)*dline_1(:,1)**6 + &
nu_aniso_hyper3(2)*dline_1(:,2)**6 + &
nu_aniso_hyper3(3)*dline_1(:,3)**6)/ dxyz_6
!
endif
!
! viscous force: mu/rho*d6uj/dx6
!
if (lvisc_hyper3_rho_nu_const_bulk) then
murho1=nu_hyper3*p%rho1 ! (=mu_hyper3/rho)
do i=1,3
p%fvisc(:,i)=p%fvisc(:,i)+murho1*p%del6u_bulk(:,i)
enddo
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_rho_nu_const_bulk')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+murho1
endif
!
! viscous force: nu_hyper3*(del6u+S.glnrho), where S_ij=d^5 u_i/dx_j^5
!
if (lvisc_hyper3_nu_const) then
p%fvisc=p%fvisc+nu_hyper3*(p%del6u+p%uij5glnrho)
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity', 'viscous heating term '// &
'is not implemented for lvisc_hyper3_nu_const')
endif
endif
if (lfirst .and. ldt) p%diffus_total3=p%diffus_total3+nu_hyper3
endif
!
! Smagorinsky viscosity
! viscous force: nu_smag*(del2u+graddivu/3+2S.glnrho)+2*gnu_smag.S
! where nu_smag=(C_smag*dxmax)**2*sqrt(2*SS^2)
!
if (lvisc_smag) then
!
if (ldensity) then
if (lnusmag_as_aux) then
p%nu_smag=f(l1:l2,m,n,inusmag)
!
! Compute gradient of p%nu_smag from f-array.
!
call grad(f,inusmag,gradnu)
else
!
! Compute nu_smag and put into a pencil
!
p%nu_smag=(C_smag*dxmax)**2.*sqrt(2.*p%sij2)
!
! Enhance nu_smag in proportion to the Mach number to power 2*nu_smag_Ma2_power,
! see e.g. Chan & Sofia (1996), ApJ, 466, 372, for a similar approach
!
if (lvisc_smag_Ma) p%nu_smag=p%nu_smag*(1.+p%Ma2**nu_smag_Ma2_power)
!
! Apply quenching term if requested
!
if (gamma_smag/=0.) p%nu_smag=p%nu_smag/sqrt(1.+gamma_smag*p%sij2)
!
endif
!
! Calculate viscous force and heating
!
call multsv_mn(p%nu_smag,p%del2u+1./3.*p%graddivu+2.*p%sglnrho,tmp)
p%fvisc=p%fvisc+tmp
if (lnusmag_as_aux) then
call multmv(p%sij,gradnu,sgradnu)
p%fvisc=p%fvisc+2.*sgradnu
endif
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*p%nu_smag*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+p%nu_smag
endif
endif
!
! Simplified Smagorinsky model (neglects the gradient of nu_smag).
! viscous force: nu_smag*(del2u+graddivu/3+2S.glnrho)
! where nu_smag=(C_smag*dxmax)**2*sqrt(2*SS^2)
!
if (lvisc_smag_simplified) then
!
if (ldensity) then
!
! Find nu_smag
!
p%nu_smag=(C_smag*dxmax)**2.*sqrt(2*p%sij2)
!
! with quenching term
!
if (gamma_smag/=0.) p%nu_smag=p%nu_smag/sqrt(1.+gamma_smag*p%sij2)
!
! Calculate viscous force
!
call multsv_mn(p%nu_smag,p%sglnrho,tmp2)
call multsv_mn(p%nu_smag,p%del2u+1./3.*p%graddivu,tmp)
p%fvisc=p%fvisc+2*tmp2+tmp
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+2*p%nu_smag*p%sij2
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+p%nu_smag
endif
endif
!
! viscous force: nu_smag*(del2u+graddivu/3+2S.glnrho)
! where nu_smag=(C_smag*dxmax)**2*sqrt(S:J)
!
if (lvisc_smag_cross_simplified) then
if (ldensity) then
p%nu_smag=(C_smag*dxmax)**2.*p%ss12
!
! Calculate viscous force
!
call multsv_mn(p%nu_smag,p%sglnrho,tmp2)
call multsv_mn(p%nu_smag,p%del2u+1./3.*p%graddivu,tmp)
p%fvisc=p%fvisc+2*tmp2+tmp
if (lpencil(i_visc_heat)) then ! Heating term not implemented
if (headtt) then
call warning('calc_pencils_viscosity','viscous heating term '//&
'is not implemented for lvisc_smag_cross_simplified')
endif
endif
if (lfirst .and. ldt) p%diffus_total=p%diffus_total+p%nu_smag
endif
endif
!
! Calculate viscouse force for slope limited diffusion
! following Rempel (2014). Here the divergence of the flux is used.
!
if (lvisc_slope_limited .and. llast) then
!
! tmp3 : divergence of flux
! tmp4 : heating, switch to turn on heating calc.
!
!p%char_speed_slope=f(l1:l2,m,n,iFF_char_c) ! old implementation
!p%fvisc=p%fvisc-f(l1:l2,m,n,iFF_div_uu:iFF_div_uu+2) ! old implementation
!
! use sld only in last sub timestep
!
if (lviscosity_heat) then
if (lcylindrical_coords .or. lspherical_coords) then
do i=1,3
call calc_slope_diff_flux(f,iux+(i-1),p,h_sld_visc,nlf_sld_visc,tmp2(:,i),div_sld_visc, &
HEAT=tmp4,HEAT_TYPE='viscose', &
FLUX1=d_sld_flux(:,1,i),FLUX2=d_sld_flux(:,2,i),FLUX3=d_sld_flux(:,3,i))
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+max(0.0,tmp4)/p%rho
enddo
if (lcylindrical_coords) then
p%fvisc(:,1)=p%fvisc(:,1)+tmp2(:,1)-(d_sld_flux(:,2,2))/x(l1:l2)
p%fvisc(:,2)=p%fvisc(:,2)+tmp2(:,2)+(d_sld_flux(:,2,1))/x(l1:l2)
p%fvisc(:,3)=p%fvisc(:,3)+tmp2(:,3)
elseif(lspherical_coords) then
if (lsld_notensor) then
p%fvisc(:,1)=p%fvisc(:,1)+tmp2(:,1)
p%fvisc(:,2)=p%fvisc(:,2)+tmp2(:,2)
p%fvisc(:,3)=p%fvisc(:,3)+tmp2(:,3)
else
p%fvisc(:,1)=p%fvisc(:,1)+tmp2(:,1)-(d_sld_flux(:,2,2)+d_sld_flux(:,3,3))/x(l1:l2)
p%fvisc(:,2)=p%fvisc(:,2)+tmp2(:,2)+(d_sld_flux(:,2,1)-d_sld_flux(:,3,3)*cotth(m))/x(l1:l2)
p%fvisc(:,3)=p%fvisc(:,3)+tmp2(:,3)+(d_sld_flux(:,3,1)+d_sld_flux(:,3,2)*cotth(m))/x(l1:l2)
endif
endif
else
do i=1,3
call calc_slope_diff_flux(f,iux+(i-1),p,h_sld_visc,nlf_sld_visc,tmp3,div_sld_visc,&
HEAT=tmp4,HEAT_TYPE='viscose')
p%fvisc(:,i)=p%fvisc(:,i) + tmp3
if (lpencil(i_visc_heat)) p%visc_heat=p%visc_heat+max(0.0,tmp4)/p%rho
enddo
endif
else
if (lcylindrical_coords .or. lspherical_coords) then
do i=1,3
call calc_slope_diff_flux(f,iux+(i-1),p,h_sld_visc,nlf_sld_visc,tmp2(:,i),div_sld_visc, &
FLUX1=d_sld_flux(:,1,i),FLUX2=d_sld_flux(:,2,i),FLUX3=d_sld_flux(:,3,i))
enddo
if (lcylindrical_coords) then
p%fvisc(:,1)=p%fvisc(:,1)+tmp2(:,1)-(d_sld_flux(:,2,2))/x(l1:l2)
p%fvisc(:,2)=p%fvisc(:,2)+tmp2(:,2)+(d_sld_flux(:,2,1))/x(l1:l2)
p%fvisc(:,3)=p%fvisc(:,3)+tmp2(:,3)
elseif(lspherical_coords) then
p%fvisc(:,1)=p%fvisc(:,1)+tmp2(:,1)-(d_sld_flux(:,2,2)+d_sld_flux(:,3,3))/x(l1:l2)
p%fvisc(:,2)=p%fvisc(:,2)+tmp2(:,2)+(d_sld_flux(:,2,1)-d_sld_flux(:,3,3)*cotth(m))/x(l1:l2)
p%fvisc(:,3)=p%fvisc(:,3)+tmp2(:,3)+(d_sld_flux(:,3,1)+d_sld_flux(:,3,2)*cotth(m))/x(l1:l2)
endif
else
do i=1,3
call calc_slope_diff_flux(f,iux+(i-1),p,h_sld_visc,nlf_sld_visc,tmp3,div_sld_visc)
p%fvisc(:,i)=p%fvisc(:,i) + tmp3
enddo
endif
endif
endif
!
! Calculate Lambda effect
!
if (llambda_effect) then
call calc_lambda(p,lambda_phi)
p%fvisc(:,iuz)=p%fvisc(:,iuz) + lambda_phi
endif
!
! Store viscous heating rate in auxiliary variable if requested.
! Just necessary immediately before writing snapshots, but how would we
! know we are?
!
if (lvisc_heat_as_aux) f(l1:l2,m,n,ivisc_heat) = p%visc_heat
!
! Store viscous force (acceleration) in auxiliary variable if requested.
!
if (lvisc_forc_as_aux) then
do i=0,2
f(l1:l2,m,n,ivisc_forc+i) = p%fvisc(:,i+1)
enddo
endif
!
! Do diagnostics related to viscosity.
!
if (ldiagnos) then
if (idiag_nusmagm/=0) call sum_mn_name(p%nu_smag,idiag_nusmagm)
if (idiag_nu_LES/=0) call sum_mn_name(p%nu_smag,idiag_nu_LES)
if (idiag_nusmagmin/=0) call max_mn_name(-p%nu_smag,idiag_nusmagmin,lneg=.true.)
if (idiag_nusmagmax/=0) call max_mn_name(p%nu_smag,idiag_nusmagmax)
if (idiag_num/=0) call sum_mn_name(p%nu,idiag_num)
if (idiag_qfviscm/=0) then
call dot(p%curlo,p%fvisc,qfvisc)
call sum_mn_name(qfvisc,idiag_qfviscm)
endif
endif
!
endsubroutine calc_pencils_viscosity
!***********************************************************************
subroutine viscosity_after_boundary(f)
!
! 19-dec-16/MR: fixed bug: m,n must be from Cdata. Added precalculation of nu_smag in f array.
! Rewritten viscous heating from slope-limited diffusion.
! 18-may-17/MR: corrected wrong order of loops for viscous heating by slope-limited diffusion.
! 16-jun-18/axel: reverted "severe bug..." of r74487; see "!AB_REM:" and "!AB_ADD"
! 16-aug-18/MR: replaced Joerns implementation of linear interpolation from staggered to normal
! grid as it was not bi-/trilinear in 2D/3D;
! added 3rd order interpolation (default)
! 03-apr-20/joern: restructured and fixed slope-limited diffusion, SLD commented out here.
!
use Sub, only: div, calc_all_diff_fluxes, grad, dot_mn, calc_sij2, &
stagger_to_base_interp_1st, stagger_to_base_interp_3rd
use General, only: reduce_grad_dim,notanumber
use DensityMethods, only: getrho
use SharedVariables, only: get_shared_variable
use Boundcond, only: update_ghosts
real, dimension(mx,my,mz,mfarray) :: f
real, dimension(nx) :: rho, tmp
logical, pointer :: lshear_rateofstrain
if (lnusmag_as_aux) then
!
! Compute nu_smag and put into tmp
!
call get_shared_variable('lshear_rateofstrain',lshear_rateofstrain,caller='viscosity_after_boundary')
do n=n1,n2; do m=m1,m2
!
! sij2 -> rho
!
call calc_sij2(f,rho,lshear_rateofstrain)
tmp=(C_smag*dxmax)**2.*sqrt(2.*rho)
!
! Enhance nu_smag in proportion to the Mach number to power 2*nu_smag_Ma2_power,
! see e.g. Chan & Sofia (1996), ApJ, 466, 372, for a similar approach
!
!!if (lvisc_smag_Ma) tmp=tmp*(1.+p%Ma2**nu_smag_Ma2_power)
!
! Apply quenching term if requested
!
if (gamma_smag/=0.) tmp=tmp/sqrt(1.+gamma_smag*rho)
!
! Put nu_smag into the f-array.
!
f(l1:l2,m,n,inusmag)=tmp
!
enddo; enddo
!
call update_ghosts(f,inusmag)
!
endif
!
! Slope limited diffusion following Rempel (2014).
! First calculating the flux in a subroutine below
! using a slope limiting procedure then storing in the
! auxiliary variables in the f array (done above).
!
! if (lvisc_slope_limited.and.lfirst) then
!
! if (lviscosity_heat) f(:,:,:,iFF_heat)=0.
!
! do j=1,3
!
! call calc_all_diff_fluxes(f,iuu+j-1,islope_limiter,h_slope_limited)
!
! do n=n1,n2; do m=m1,m2
!
! Divergence of flux = force/mass = acceleration.
!
! call div(f,iFF_diff,f(l1:l2,m,n,iFF_div_uu+j-1),.true.)
!
! Heating term.
!
! if (lviscosity_heat) then
!
! grad(rho) or grad(lnrho) (reference state?)
!
! call grad(f,ilnrho,gr)
!
! compactify the non-zero components in the first dimensionality elements of gr
!
! call reduce_grad_dim(gr)
! if (ldensity_nolog) then
! call getrho(f(:,m,n,irho),rho)
! do k=1,dimensionality
!
! now grad(ln(rho))
!
! gr(:,k)=gr(:,k)/rho
! enddo
! endif
!
! grad(u_j)
!
! call grad(f,iuu+j-1,guj)
!
! compactify the non-zero components in the first dimensionality elements of guj
!
! call reduce_grad_dim(guj)
!
! grad(ln(rho))*u_j+grad(u_j)=grad(rho*u_j)/rho
!
! do k=1,dimensionality
! guj(:,k)=gr(:,k)*f(l1:l2,m,n,iuu+j-1)+guj(:,k)
! enddo
!
! loop inside has length dimensionality
! \partial_k(rho*u_j) f_jk/rho = energy/time/mass (summation over j by loop)
!
!
! do k=1,dimensionality
! if (dim_mask(k) == 1) &
! f(l1:l2,m,n,iFF_heat)=f(l1:l2,m,n,iFF_heat)+guj(:,k)* &
! (f(l1:l2,m,n,iFF_diff1-1+k)+f(l1+1:l2+1,m,n,iFF_diff1-1+k))/2.
! if (dim_mask(k) == 2) &
! f(l1:l2,m,n,iFF_heat)=f(l1:l2,m,n,iFF_heat)+guj(:,k)* &
! (f(l1:l2,m,n,iFF_diff1-1+k)+f(l1:l2,m+1,n,iFF_diff1-1+k))/2.
! if (dim_mask(k) == 3) &
! f(l1:l2,m,n,iFF_heat)=f(l1:l2,m,n,iFF_heat)+guj(:,k)* &
! (f(l1:l2,m,n,iFF_diff1-1+k)+f(l1:l2,m,n+1,iFF_diff1-1+k))/2.
! call stagger_to_base_interp_3rd(f(:,:,:,iFF_diff1-1+k),m,n,tmp)
!f(:,m ,n,iFF_diff1-1+k)=n+(/1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12./)
! call stagger_to_base_interp_1st(f(:,:,:,iFF_diff1-1+k),m,n,tmp)
! f(l1:l2,m,n,iFF_heat)=f(l1:l2,m,n,iFF_heat)+guj(:,k)*tmp
!print*, 'k,tmp=', k, tmp
! enddo
! call dot_mn(guj(:,1:dimensionality),f(l1:l2,m,n,iFF_diff1:iFF_diff2), &
! f(l1:l2,m,n,iFF_heat),ladd=.true.)
!
! no cooling admitted (Why?)
!
! f(l1:l2,m,n,iFF_heat)=min(f(l1:l2,m,n,iFF_heat),0.)
! endif
! enddo; enddo
! enddo
!maxh=maxval(abs(f(l1:l2,m1:m2,n1:n2,iFF_heat)))
!if (ldiagnos.and.maxh>1.) print*, 'heat:', iproc, maxh !, minval(f(l1:l2,m1:m2,n1:n2,iFF_heat))
!if (ldiagnos) print*, 'grhouj:', iproc, maxval(guj(:,1:dimensionality)), minval(guj(:,1:dimensionality))
!f(l1:l2,m1:m2,n1:n2,iFF_heat)=0. !!!
! endif
endsubroutine viscosity_after_boundary
!***********************************************************************
subroutine getnu_non_newtonian(gdotsqr,nu_effective,gradnu_effective)
!
! Outputs non-newtonian viscosity for an input strain-rate squared
!
use Sub, only : step,der_step
real, dimension(nx) :: gdotsqr
real,dimension(nx) :: nu_effective,gradnu_effective
!
select case(nnewton_type)
case('carreau')
nu_effective= nu_infinity+ &
(nu0-nu_infinity)/(1+(non_newton_lambda**2)*gdotsqr)**((1.-carreau_exponent)/2.)
gradnu_effective= (non_newton_lambda**2)*0.5*(carreau_exponent-1.)* &
(nu0-nu_infinity)/(1+(non_newton_lambda**2)*gdotsqr)**(1.5-carreau_exponent/2.)
case('step')
nu_effective=nu0+(nu_infinity-nu0)*step(gdotsqr,nnewton_tscale**2,nnewton_step_width)
gradnu_effective=(nu_infinity-nu0)*der_step(gdotsqr,nnewton_tscale**2,nnewton_step_width)
case default
write(*,*) 'nnewton_type=',nnewton_type
call fatal_error('viscosity,getnu_non_newtonian:','select nnewton_type')
endselect
!
endsubroutine getnu_non_newtonian
!***********************************************************************
subroutine calc_viscosity(f)
!
real, dimension(mx,my,mz,mfarray), intent(in) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine calc_viscosity
!***********************************************************************
subroutine get_gradnu(tmp,ljumpx,ljumpr,p,gradnu)
!
! Calculate grad nu for vicosity jump in different
! coordinate systems
!
! 20-jun-08/wlad :: coded
!
real, dimension(nx) ,intent(in) :: tmp
real, dimension(nx,3),intent(out) :: gradnu
logical, intent (in) :: ljumpx,ljumpr
type (pencil_case) :: p
!
if (ljumpx.or. &
(ljumpr.and.(lcylindrical_coords.or. &
lspherical_coords ))) then
gradnu(:,1)=tmp ; gradnu(:,2)=0 ; gradnu(:,3)=0
elseif (ljumpr.and.lcylinder_in_a_box) then
gradnu(:,1)=tmp*x(l1:l2)*p%rcyl_mn1
gradnu(:,2)=tmp*y( m )*p%rcyl_mn1
elseif (ljumpr.and.lsphere_in_a_box) then
gradnu(:,1)=tmp*x(l1:l2)*p%r_mn1
gradnu(:,2)=tmp*y( m )*p%r_mn1
gradnu(:,3)=tmp*z( n )*p%r_mn1
endif
!
endsubroutine get_gradnu
!***********************************************************************
subroutine calc_viscous_heat(df,p,Hmax)
!
! Add viscous heating term to the right hand side of entropy equation.
!
! 20-nov-02/tony: coded
!
use Sub, only: cubic_step
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
real, dimension (nx),intent(inout) :: Hmax
!
! Add viscous heat (which has units of energy/mass) to the RHS
! of the entropy (both with and without pretend_lnTT), or of
! the temperature equation. Divide by cv if pretend_lnTT.
!
if (lentropy) then
if (pretend_lnTT) then
df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + p%cv1*p%TT1*p%visc_heat
else
df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + p%TT1*p%visc_heat
endif
else if (ltemperature) then
if (ltemperature_nolog) then
df(l1:l2,m,n,iTT) = df(l1:l2,m,n,iTT) + p%cv1*p%visc_heat
else
if (lno_visc_heat_zbound) then
df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + &
p%cv1*p%TT1*p%visc_heat*cubic_step(z(n),no_visc_heat_z0,no_visc_heat_zwidth)
else
df(l1:l2,m,n,ilnTT) = df(l1:l2,m,n,ilnTT) + p%cv1*p%TT1*p%visc_heat
endif
endif
else if (lthermal_energy) then
if (lstratz) then
df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + p%rho * p%visc_heat / eth0z(n)
else
df(l1:l2,m,n,ieth) = df(l1:l2,m,n,ieth) + p%rho * p%visc_heat
endif
endif
!
! Calculate maximum heating (for time step constraint), so it is
! only done on the first of the 3 substeps.
!
if (lfirst .and. ldt) Hmax=Hmax+p%visc_heat
!
endsubroutine calc_viscous_heat
!***********************************************************************
subroutine calc_viscous_force(df,p)
!
! Calculate viscous force term for right hand side of equation of motion.
!
! 20-nov-02/tony: coded
! 9-jul-04/nils: added Smagorinsky viscosity
!
use Diagnostics, only: sum_mn_name, max_mn_name, xysum_mn_name_z, &
yzsum_mn_name_x, zsum_mn_name_xy, max_mn_name, phisum_mn_name_rz, &
integrate_mn_name
use Sub, only: cross, dot2, mult_mat_vv
!
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
intent (in) :: p
intent (inout) :: df
!
real, dimension (nx) :: Reshock, fvisc2, uus, tmp
real, dimension (nx,3):: nuD2uxb, fluxv
real, dimension (nx) :: diffus_nu, diffus_nu2, diffus_nu3
integer :: i
!
! Add viscosity to equation of motion
!
if (lmagfield_nu) then
do i=1,3
df(l1:l2,m,n,iux+i-1) = df(l1:l2,m,n,iux+i-1) + p%fvisc(:,i)/(1.+p%b2/meanfield_nuB**2)
enddo
else
df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) + p%fvisc
endif
!
! Calculate max total diffusion coefficient for timestep calculation etc.
!
if (lfirst.and.ldt) then
diffus_nu =p%diffus_total *dxyz_2
diffus_nu2=p%diffus_total2*dxyz_4
if (ldynamical_diffusion .and. lvisc_hyper3_mesh) then
diffus_nu3 = p%diffus_total3 * sum(abs(dline_1),2)
else
diffus_nu3=p%diffus_total3*dxyz_6
endif
maxdiffus=max(maxdiffus,diffus_nu)
maxdiffus2=max(maxdiffus2,diffus_nu2)
maxdiffus3=max(maxdiffus3,diffus_nu3)
if (ldiagnos.and.idiag_dtnu/=0) &
call max_mn_name(diffus_nu/cdtv,idiag_dtnu,l_dt=.true.)
if (ldiagnos.and.idiag_dtnu3/=0) &
call max_mn_name(diffus_nu3/cdtv3,idiag_dtnu3,l_dt=.true.)
endif
!
! Diagnostic output
!
if (ldiagnos) then
if (idiag_nu_tdep/=0) call sum_mn_name(spread(nu_tdep,1,nx),idiag_nu_tdep)
! if (lvisc_smag_simplified) then
! if (ldensity) nu_smag=(C_smag*dxmax)**2.*sqrt(2*p%sij2)
! endif
! if (lvisc_smag_cross_simplified) then
! if (ldensity) nu_smag=(C_smag*dxmax)**2.*p%ss12
! endif
if (idiag_meshRemax/=0) call max_mn_name(sqrt(p%u2(:))*dxmax_pencil/p%diffus_total,idiag_meshRemax)
if (idiag_Reshock/=0) then
Reshock(:) = 0.
where (abs(p%shock) > tini)
Reshock = dxmax_pencil*sqrt(p%u2)/(nu_shock*p%shock)
endwhere
call max_mn_name(Reshock,idiag_Reshock)
endif
if (idiag_Sij2m/=0) call sum_mn_name(p%sij2,idiag_Sij2m)
!
! Viscous heating for Smagorinsky viscosity.
!
if (idiag_epsK_LES/=0) then
if (lvisc_smag_simplified) then
call sum_mn_name(2*p%nu_smag*p%rho*p%sij2,idiag_epsK_LES)
else if (lvisc_smag_cross_simplified) then
call sum_mn_name(2*p%nu_smag*p%rho*p%sij2,idiag_epsK_LES)
endif
endif
!
! Calculate sijoiojm = S_{ij} o_i o_j.
!
if (idiag_sijoiojm/=0) then
call mult_mat_vv(p%sij,p%oo,p%oo,tmp)
call sum_mn_name(tmp,idiag_sijoiojm)
endif
!
! correlation of viscous term with b-field; relevant for MTA
! (MTA: minimal tau approximation)
!
if (lmagnetic) then
if (idiag_nuD2uxbxm/=0.or.idiag_nuD2uxbym/=0.or.idiag_nuD2uxbzm/=0) then
call cross(p%fvisc,p%bb,nuD2uxb)
call sum_mn_name(nuD2uxb(:,1),idiag_nuD2uxbxm)
call sum_mn_name(nuD2uxb(:,2),idiag_nuD2uxbym)
call sum_mn_name(nuD2uxb(:,3),idiag_nuD2uxbzm)
endif
endif
endif
!
! For slope-limted diffusion viscocity diagnostics need to be written
! out at the last time step
!
if (ldiagnos) then
if (idiag_fviscm/=0 .or. idiag_fviscmax/=0 .or. idiag_fviscrmsx/=0) &
call dot2(p%fvisc,fvisc2)
if (idiag_fviscm/=0) call sum_mn_name(fvisc2,idiag_fviscm,lsqrt=.true.)
if (idiag_ufviscm/=0) &
call sum_mn_name(p%uu(:,1)*p%fvisc(:,1)+ &
p%uu(:,2)*p%fvisc(:,2)+ &
p%uu(:,3)*p%fvisc(:,3),idiag_ufviscm)
if (idiag_fviscmin/=0) call max_mn_name(-p%fvisc,idiag_fviscmin,lneg=.true.)
if (idiag_fviscmax/=0) call max_mn_name(fvisc2,idiag_fviscmax,lsqrt=.true.)
if (idiag_fviscrmsx/=0) call sum_mn_name(xmask_vis*fvisc2,idiag_fviscrmsx,lsqrt=.true.)
if (idiag_visc_heatm/=0) call sum_mn_name(p%visc_heat,idiag_visc_heatm)
if (idiag_epsK/=0) call sum_mn_name(p%visc_heat*p%rho,idiag_epsK)
if (idiag_epsKint/=0) call integrate_mn_name(p%visc_heat*p%rho,idiag_epsKint)
endif
!
! 1D-averages.
!
if (l1davgfirst) then
if (idiag_fviscmz/=0) &
call xysum_mn_name_z(-2.*p%rho*nu*( &
p%uu(:,1)*p%sij(:,1,3)+ &
p%uu(:,2)*p%sij(:,2,3)+ &
p%uu(:,3)*p%sij(:,3,3)),idiag_fviscmz)
if (idiag_fviscsmmz/=0) &
call xysum_mn_name_z(-2.*p%rho*p%nu_smag*( &
p%uu(:,1)*p%sij(:,1,3)+ &
p%uu(:,2)*p%sij(:,2,3)+ &
p%uu(:,3)*p%sij(:,3,3)),idiag_fviscsmmz)
if (idiag_epsKmz/=0) &
call xysum_mn_name_z(p%visc_heat*p%rho,idiag_epsKmz)
if (idiag_fviscmx/=0) &
call yzsum_mn_name_x(-2.*p%rho*nu*( &
p%uu(:,1)*p%sij(:,1,1)+ &
p%uu(:,2)*p%sij(:,2,1)+ &
p%uu(:,3)*p%sij(:,3,1)),idiag_fviscmx)
if (idiag_numx/=0) &
call yzsum_mn_name_x(p%nu,idiag_numx)
if (idiag_viscforcezmz/=0) &
call xysum_mn_name_z(p%rho*p%fvisc(:,3),idiag_viscforcezmz)
if (idiag_viscforcezupmz/=0) then
where (p%uu(:,3) > 0.)
uus = p%rho*p%fvisc(:,3)
elsewhere
uus=0.
endwhere
call xysum_mn_name_z(uus,idiag_viscforcezupmz)
endif
if (idiag_viscforcezdownmz/=0) then
where (p%uu(:,3) < 0.)
uus = p%rho*p%fvisc(:,3)
elsewhere
uus=0.
endwhere
call xysum_mn_name_z(uus,idiag_viscforcezdownmz)
endif
endif
!
! 2D-averages.
!
if (l2davgfirst) then
if (idiag_fviscmxy/=0) then
if (lvisc_sqrtrho_nu_const) then
call zsum_mn_name_xy(-2.*sqrt(p%rho)*nu*( &
p%uu(:,1)*p%sij(:,1,1)+p%uu(:,2)*p%sij(:,2,1)+ &
p%uu(:,3)*p%sij(:,3,1)),idiag_fviscmxy)
elseif (lvisc_rho_nu_const) then
call zsum_mn_name_xy(-2.*nu*( &
p%uu(:,1)*p%sij(:,1,1)+p%uu(:,2)*p%sij(:,2,1)+ &
p%uu(:,3)*p%sij(:,3,1)),idiag_fviscmxy)
elseif (lvisc_nu_profy_bound) then
call zsum_mn_name_xy(-2.*p%rho*pnu*( &
p%uu(:,1)*p%sij(:,1,1)+p%uu(:,2)*p%sij(:,2,1)+ &
p%uu(:,3)*p%sij(:,3,1)),idiag_fviscmxy)
else
call zsum_mn_name_xy(-2.*p%rho*nu*( &
p%uu(:,1)*p%sij(:,1,1)+p%uu(:,2)*p%sij(:,2,1)+ &
p%uu(:,3)*p%sij(:,3,1)),idiag_fviscmxy)
endif
endif
if (idiag_fviscsmmxy/=0) call zsum_mn_name_xy(-2.*p%rho*p%nu_smag*( &
p%uu(:,1)*p%sij(:,1,1)+p%uu(:,2)*p%sij(:,2,1)+ &
p%uu(:,3)*p%sij(:,3,1)),idiag_fviscsmmxy)
if (idiag_fviscymxy/=0) then
if (lyang) then
fluxv(:,1)=0.
fluxv(:,2)=p%uu(:,1)*p%sij(:,1,2)+p%uu(:,2)*p%sij(:,2,2)+p%uu(:,3)*p%sij(:,3,2)
fluxv(:,3)=p%uu(:,1)*p%sij(:,1,3)+p%uu(:,2)*p%sij(:,2,3)+p%uu(:,3)*p%sij(:,3,3)
call zsum_mn_name_xy(fluxv,idiag_fviscymxy,(/0,1,0/),-2.*p%rho*nu)
else
call zsum_mn_name_xy(-2.*p%rho*nu*( &
p%uu(:,1)*p%sij(:,1,2)+p%uu(:,2)*p%sij(:,2,2)+ &
p%uu(:,3)*p%sij(:,3,2)),idiag_fviscymxy)
endif
endif
if (idiag_fviscrsphmphi/=0) then
fluxv(:,1)=p%uu(:,1)*p%sij(:,1,1)+p%uu(:,2)*p%sij(:,2,1)+p%uu(:,3)*p%sij(:,3,1)
fluxv(:,2)=p%uu(:,1)*p%sij(:,1,2)+p%uu(:,2)*p%sij(:,2,2)+p%uu(:,3)*p%sij(:,3,2)
fluxv(:,3)=p%uu(:,1)*p%sij(:,1,3)+p%uu(:,2)*p%sij(:,2,3)+p%uu(:,3)*p%sij(:,3,3)
call phisum_mn_name_rz(-2.*p%rho*nu*(fluxv(:,1)*p%evr(:,1)+ &
fluxv(:,2)*p%evr(:,2)+fluxv(:,3)*p%evr(:,3)),idiag_fviscrsphmphi)
endif
endif
!
endsubroutine calc_viscous_force
!***********************************************************************
subroutine calc_visc_heat_ppd(df,p)
!
! Calculates viscous dissipation term from D'Angelo et al. (2003)
!
! 03/08 steveb
!
use Sub, only: get_radial_distance
use Gravity, only: acceleration
!
real, dimension (mx,my,mz,mvar) :: df
type (pencil_case) :: p
!
real, dimension (nx) :: diss
real, dimension (nx) :: rr_sph,rr_cyl,g_r,OO2
!
! 'Dissipative' heating term Y=9/4 \Sigma \nu \Omega_K^2
! Need to get correct angular velocity first
!
call get_radial_distance(rr_sph,rr_cyl)
call acceleration(g_r)
! should generalize this to non-cartesian; see hydro L1937
OO2=max(-g_r/rr_cyl,0.)
!
! only works for 2-D r-phi disks
! this goes into entropy equation, which divides by \rho, so no density here
! this is really surface density, though
!
if (nzgrid==1) then
diss = 9./4.*nu*OO2
df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + p%TT1*diss
!
else
call fatal_error("calc_visc_heat_ppd", &
"dissipation only implemented for 2d-disk")
endif
!
endsubroutine calc_visc_heat_ppd
!***********************************************************************
subroutine getnu(p,nu_input,nu_pencil,ivis)
!
! Return the viscosity (by value)
!
! 14-aug-08/kapelrud: coded
!
type (pencil_case), optional :: p
real, optional, intent(out) :: nu_input
real, dimension(nx), optional, intent(out) :: nu_pencil
character (len=labellen), optional :: ivis
!
if (present(nu_input)) nu_input=nu
if (present(ivis)) ivis=ivisc(1)
if (present(nu_pencil)) then
if (.not. present(p)) then
call fatal_error('getnu',&
'p must be present in call to getnu when nu_pencil is present!')
endif
if (lvisc_simplified) then
nu_pencil=nu
elseif(lvisc_mixture) then
nu_pencil=nu
elseif (lvisc_rho_nu_const) then
nu_pencil=nu*p%rho1
elseif (lvisc_rho_nu_const_bulk) then
nu_pencil=zeta*p%rho1
elseif (lvisc_nu_cspeed) then
nu_pencil=nu*exp(0.5*p%lnTT)
elseif (lvisc_mu_cspeed) then
nu_pencil=nu*exp(0.5*p%lnTT)*p%rho1
elseif (lvisc_nu_const) then
nu_pencil=nu
else
call fatal_error('getnu','No such ivisc implemented!')
endif
endif
!
endsubroutine getnu
!***********************************************************************
subroutine dynamical_viscosity(uc)
!
! Dynamically set viscosity coefficient given fixed mesh Reynolds number.
!
! 27-jul-11/ccyang: coded
!
! Input Argument
! uc
! Characteristic velocity of the system.
!
real, intent(in) :: uc
!
! Hyper-viscosity coefficient
!
if (nu_hyper3 /= 0.0) nu_hyper3 = pi5_1 * uc * dxmax**5 / re_mesh
if (nu_hyper3_mesh /= 0.0) nu_hyper3_mesh = pi5_1 * uc / re_mesh / sqrt(real(dimensionality))
!
endsubroutine dynamical_viscosity
!***********************************************************************
subroutine split_update_viscosity(f)
!
! Update the velocity by integrating the operator split viscous terms.
!
! 22-aug-13/ccyang: coded.
!
use ImplicitDiffusion, only: integrate_diffusion
!
real, dimension(mx,my,mz,mfarray), intent(inout) :: f
!
if (limplicit_viscosity) &
call integrate_diffusion(get_viscosity_implicit, f, iux, iuz)
!
endsubroutine split_update_viscosity
!***********************************************************************
subroutine get_viscosity_implicit(ndc, diffus_coeff, iz)
!
! Gets the diffusion coefficient along a given pencil for the implicit algorithm.
!
! 22-aug-13/ccyang: coded.
!
integer, intent(in) :: ndc
real, dimension(ndc), intent(out) :: diffus_coeff
integer, intent(in), optional :: iz
!
if (lvisc_simplified) then
diffus_coeff = nu
else
diffus_coeff = 0.0
endif
!
endsubroutine get_viscosity_implicit
!***********************************************************************
subroutine calc_lambda(p,div_lambda)
!
! Calculates the lambda effect
!
! 20-apr-10/dhruba: coded
! If lKit_Olem is true the lambda coefficients depends on dsdr which must be
! incorporated below.
!
!
use cdata, only: Omega
!
real,dimension(nx) :: div_lambda,lomega,dlomega_dr,dlomega_dtheta, &
lver,lhor,dlver_dr,dlhor_dtheta
type (pencil_case) :: p
!
lomega=p%uu(:,3)/(sinth(m)*x(l1:l2))+Omega
!
dlomega_dr=(x(l1:l2)*p%uij(:,3,1)-p%uu(:,3))/ &
(sinth(m)*x(l1:l2)*x(l1:l2))
dlomega_dtheta=(p%uij(:,3,2)*x(l1:l2)-p%uu(:,3)*cotth(m))/ &
(sinth(m)*x(l1:l2)*x(l1:l2))
!
lver = -(Lambda_V0*LV0_rprof(l1:l2)+Lambda_V1*sinth(m)*sinth(m) &
*LV1_rprof(l1:l2) )
lhor = -Lambda_H1*sinth(m)*sinth(m)*LH1_rprof(l1:l2)
!
dlver_dr = -(Lambda_V0*der_LV0_rprof(l1:l2)+Lambda_V1 &
*sinth(m)*sinth(m)*der_LV1_rprof(l1:l2))
dlhor_dtheta = -Lambda_H1*LH1_rprof(l1:l2)*2.*costh(m)*sinth(m)/x(l1:l2)
!
div_lambda = lver*(sinth(m)*lomega*p%glnrho(:,1) &
+3.*sinth(m)*lomega/x(l1:l2) + sinth(m)*dlomega_dr) &
+lomega*sinth(m)*dlver_dr + lhor*(costh(m)*lomega*p%glnrho(:,2) &
-sinth(m)*lomega/x(l1:l2) + 2.*cotth(m)*costh(m)*lomega/x(l1:l2) &
+costh(m)*dlomega_dtheta) + lomega*costh(m)*dlhor_dtheta
!
endsubroutine calc_lambda
!***********************************************************************
subroutine pushdiags2c(p_diag)
integer, parameter :: n_diags=0
integer(KIND=ikind8), dimension(:) :: p_diag
call keep_compiler_quiet(p_diag)
endsubroutine pushdiags2c
!***********************************************************************
subroutine pushpars2c(p_par)
use Syscalls, only: copy_addr
integer, parameter :: n_pars=2
integer(KIND=ikind8), dimension(n_pars) :: p_par
call copy_addr(nu,p_par(1))
call copy_addr(zeta,p_par(2))
endsubroutine pushpars2c
!***********************************************************************
!
endmodule Viscosity