! $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 ! !*************************************************************** module Viscosity use Cparam use Cdata use Messages use Density implicit none include 'viscosity.h' ! real :: nu=0. character (len=labellen) :: ivisc='nu-const' real :: nu_var, q_DJO=2., t0_DJO=0., nuf_DJO,ti_DJO=1.,tf_DJO integer :: dummy integer :: idiag_dtnu=0 ! dummy logical logical :: lvisc_first=.false. ! run parameters namelist /viscosity_run_pars/ ivisc,q_DJO,t0_DJO,nu,nuf_DJO,ti_DJO,tf_DJO ! Not implemented but needed for bodged implementation in hydro integer :: idiag_epsK=0 contains !*********************************************************************** subroutine register_viscosity() ! ! 19-nov-02/tony: coded ! use Cdata use Mpicomm use Sub ! if ((ip<=8) .and. lroot) then print*, 'register_viscosity: constant viscosity' endif ! ! Identify version number. ! if (lroot) call svn_id( & "$Id$") ! endsubroutine register_viscosity !*********************************************************************** subroutine initialize_viscosity() ! ! 20-nov-02/tony: coded ! use Cdata ! if (nu /= 0. .and. ((ivisc=='nu-const') .or. (ivisc=='nu-DJO'))) then lneed_sij=.true. endif ! endsubroutine initialize_viscosity !*********************************************************************** 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 Cdata use Diagnostics use FArrayManager, only: farray_index_append ! logical :: lreset logical, optional :: lwrite integer :: iname ! ! reset everything in case of reset ! (this needs to be consistent with what is defined above!) ! if (lreset) then idiag_dtnu=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),'dtnu',idiag_dtnu) enddo ! ! write column where which ionization variable is stored ! if (present(lwrite)) then if (lwrite) then call farray_index_append('i_dtnu',idiag_dtnu) call farray_index_append('ihyper',ihyper) call farray_index_append('ishock',ishock) call farray_index_append('itest',0) endif endif ! endsubroutine rprint_viscosity !*********************************************************************** subroutine calc_viscosity(f) ! real, dimension (mx,my,mz,mfarray) :: f ! call keep_compiler_quiet(f) ! endsubroutine calc_viscosity !*********************************************************************** subroutine calc_viscous_heat(f,df,glnrho,divu,rho1,cs2,TT1,shock) ! ! calculate viscous heating term for right hand side of entropy equation ! ! 20-nov-02/tony: coded ! use Cdata use Mpicomm use Sub real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (nx) :: rho1,TT1,cs2 real, dimension (nx) :: sij2, divu, shock real, dimension (nx,3) :: glnrho ! ! traceless strainmatrix squared ! call multm2(sij,sij2) ! select case (ivisc) case ('simplified', '0') if (headtt) print*,'no heating: ivisc=',ivisc heat=0. case ('rho_nu-const', '1') if (headtt) print*,'viscous heating: ivisc=',ivisc heat=2.*nu*sij2*rho1 case ('nu-const', '2') if (headtt) print*,'viscous heating: ivisc=',ivisc heat=2.*nu*sij2 case default if (lroot) print*,'ivisc=',trim(ivisc),' -- this could never happen' call stop_it("") !" endselect df(l1:l2,m,n,iss) = df(l1:l2,m,n,iss) + TT1*heat if (lfirst.and.ldt) Hmax=Hmax+heat endsubroutine calc_viscous_heat !*********************************************************************** subroutine calc_viscous_force(f,df,glnrho,divu,rho,rho1,shock,gshock,bij) ! use Cdata use Diagnostics use Mpicomm use Sub real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (nx,3,3) :: bij real, dimension (nx,3) :: glnrho,del2u, graddivu,fvisc,sglnrho,gshock real, dimension (nx) :: murho1,rho1,divu,shock,rho real :: b integer :: i intent (in) :: f, glnrho, rho1,rho intent (out) :: df,shock,gshock ! ! viscosity operator ! rho1 is pre-calculated in equ ! ! if (nu /= 0.) then shock=0. gshock=0. select case (ivisc) case ('nu-const') ! ! viscous force: nu*(del2u+graddivu/3+2S.glnrho) ! -- the correct expression for nu=const ! if (headtt) print*,'VISC_VAR: viscous force: nu*(del2u+graddivu/3+2S.glnrho)' call del2v_etc(f,iuu,del2u,GRADDIV=graddivu) if (ldensity) then call multmv(sij,glnrho,sglnrho) fvisc=2*nu*sglnrho+nu*(del2u+1./3.*graddivu) diffus_nu=diffus_nu+nu*dxyz_2 else if (lfirstpoint) & print*,"ldensity better be .true. for ivisc='nu-const'" endif case ('nu-DJO') ! Ditlevsen, Jensen, and Olesen (Phys Rev Lett) show that ! E=k^q*psi((t*k^(3+q)/2),nu t^{-(1-q)/(3+q)}), t_code=t-t0 ! Viscosity is choosen here such that second argument remains constant ! with time. q and t0 must be guessed a priori. q=3.67 might be a good ! choice. ! if nu is given in namelist this will correspond to nu at t=1. ! alternatively, nu_min=nu(t_max) and t_max can be set in namelist ! ! viscous force: nu*(del2u+graddivu/3+2S.glnrho) ! -- the correct expression for nu=const ! b = - (1.-q_DJO)/(3.+q_DJO) if (t.lt.ti_DJO) then nu_var = nu else nu_var= nu*((ti_DJO + t0_DJO)/(t + t0_DJO))**b endif tf_DJO = (nu/nuf_DJO)**(1./b)*(ti_DJO + t0_DJO) - t0_DJO if (lroot.and.lfirstpoint.and.lout) call write_viscosity if (headtt) print*,'Using DJO variable viscosity. with ' if (headtt.and.(ip.lt.10)) then ! print*,'VISC_VAR: nu, nuf,ti,tf,t0,q,nu_var' ! print*,nu,nuf_DJO,ti_DJO,tf_DJO,t0_DJO,q_DJO,nu_var endif if (headtt) print*,'viscous force: nu_var*(del2u+graddivu/3+2S.glnrho)' call del2v_etc(f,iuu,del2u,GRADDIV=graddivu) if (ldensity) then call multmv(sij,glnrho,sglnrho) fvisc=2.*nu_var*sglnrho+nu_var*(del2u+1./3.*graddivu) diffus_nu=diffus_nu+nu_var*dxyz_2 else if (lfirstpoint) & print*,"ldensity better be .true. for ivisc=",ivisc endif case ('nu-DJO-nobulk') b = - (1.-q_DJO)/(3.+q_DJO) if (t.lt.ti_DJO) then nu_var = nu else nu_var= nu*((ti_DJO + t0_DJO)/(t + t0_DJO))**b endif tf_DJO = (nu/nuf_DJO)**(1./b)*(ti_DJO + t0_DJO) - t0_DJO if (lroot.and.lfirstpoint.and.lout) call write_viscosity if (headtt.and.(ip.lt.6)) then print*,'VISC_VAR: nu, nuf,ti,tf,t0,q' print*,nu,nuf_DJO,ti_DJO,tf_DJO,t0_DJO,q_DJO endif if (headtt) print*,'VARIABLE viscous force: nu_var*del2v' call del2v(f,iuu,del2u) fvisc=nu_var *del2u case default ! ! Catch unknown values ! if (lroot) print*, 'Program compiled with VISCOSITY = visc_var ' if (lroot) print*, 'No such value for ivisc: ', trim(ivisc) call stop_it('calc_viscous_forcing') endselect df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)+fvisc ! ! set viscous time step ! if (ldiagnos.and.idiag_dtnu/=0) then ! Note: cdtvDim is no longer used ! call max_mn_name(spread(nu,1,nx)/dxmin**2/cdtvDim,idiag_dtnu,l_dt=.true.) call max_mn_name(spread(nu,1,nx)/dxmin**2,idiag_dtnu,l_dt=.true.) endif ! end subroutine calc_viscous_force !*********************************************************************** subroutine write_viscosity open(1,file=trim(datadir)//'/visc_var.dat',position='append') write(1,*) it,t,nu_var close(1) end subroutine !*********************************************************************** endmodule Viscosity