! $Id$ ! ! This module is the dummy for the SGS_hydro module ! in which e.g the SGS force or heat is ! computed. ! !** 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 :: lSGS_hydro = .true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 9 ! COMMUNICATED AUXILIARIES 9 ! ! PENCILS PROVIDED SGS_force(3) ! PENCILS PROVIDED SGS_heat ! !*************************************************************** module SGS_hydro ! use Cdata use General, only: keep_compiler_quiet use Messages, only: fatal_error ! implicit none ! include 'SGS_hydro.h' ! integer :: itauSGSRey=0, itauSGSMax=0, iuusmooth=0, iaasmooth=0, iSGS_heat=0, iSGS_force=0 ! real :: cReyStress=1., cMaxStress=1. logical :: lSGS_heat_as_aux=.false., lSGS_forc_as_aux=.false. namelist /SGS_hydro_run_pars/ cReyStress, cMaxStress, & lSGS_heat_as_aux, lSGS_forc_as_aux integer :: idiag_fSGSm=0, idiag_fSGSrmsx=0 real, dimension(7,7,7) :: smth_kernel contains !*********************************************************************** subroutine register_SGS_hydro ! ! 19-nov-02/tony: coded ! use Messages, only: svn_id use FArrayManager, only: farray_register_auxiliary ! ! Identify version number. ! if (lroot) call svn_id( & "$Id$") ! ! Register SGS Reynolds stress as auxilliary variable. ! call farray_register_auxiliary('uusmooth',iuusmooth,vector=3,communicated=.true.) if (lmagnetic) call farray_register_auxiliary('aasmooth',iaasmooth,vector=3,communicated=.true.) ! ! Register SGS Maxwell stress as auxilliary variable. ! call farray_register_auxiliary('tauSGSRey',itauSGSRey,vector=6,communicated=.true.) if (lmagnetic) call farray_register_auxiliary('tauSGSMax',itauSGSMax,vector=6,communicated=.true.) print*, 'tauSGSRey=', itauSGSRey print*, 'uusmooth=', iuusmooth print*, 'tauSGSMax=', itauSGSMax print*, 'aasmooth=', iaasmooth ! ! Register an extra aux slot for dissipation rate if requested (so ! SGS_heat is written to snapshots and can be easily analyzed later). ! if (lSGS_heat_as_aux) then call farray_register_auxiliary('SGS_heat',iSGS_heat) aux_var(aux_count)=',SGS_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 SGS force (accelaration) if requested (so ! visc_forc is written to snapshots and can be easily analyzed later). ! if (lSGS_forc_as_aux) then call farray_register_auxiliary('SGS_force',iSGS_force,vector=3) aux_var(aux_count)=',SGS_force' if (naux+naux_com < maux+maux_com) aux_var(aux_count)=trim(aux_var(aux_count))//' $' aux_count=aux_count+3 endif ! endsubroutine register_SGS_hydro !*********************************************************************** subroutine initialize_SGS_hydro ! use Sub, only: smoothing_kernel ! if (.not.ldensity) & call fatal_error('initialize_SGS_hydro','density needed') ! call smoothing_kernel(smth_kernel) ! ! Rectifying the boundary conditions. ! bcx12(iuusmooth:iuusmooth+2,:)=bcx12(iux:iuz,:) bcy12(iuusmooth:iuusmooth+2,:)=bcy12(iux:iuz,:) bcz12(iuusmooth:iuusmooth+2,:)=bcz12(iux:iuz,:) if (.not.lperi(1)) bcx12(itauSGSRey:itauSGSRey+5,:)='tay' if (.not.lperi(2)) bcy12(itauSGSRey:itauSGSRey+5,:)='tay' if (.not.lperi(3)) bcz12(itauSGSRey:itauSGSRey+5,:)='tay' if (lmagnetic) then bcx12(iaasmooth:iaasmooth+2,:)=bcx12(iax:iaz,:) bcy12(iaasmooth:iaasmooth+2,:)=bcy12(iax:iaz,:) bcz12(iaasmooth:iaasmooth+2,:)=bcz12(iax:iaz,:) if (.not.lperi(1)) bcx12(itauSGSMax:itauSGSMax+5,:)='tay' if (.not.lperi(2)) bcy12(itauSGSMax:itauSGSMax+5,:)='tay' if (.not.lperi(3)) bcz12(itauSGSMax:itauSGSMax+5,:)='tay' endif endsubroutine initialize_SGS_hydro !*********************************************************************** subroutine read_SGS_hydro_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=SGS_hydro_run_pars, IOSTAT=iostat) ! endsubroutine read_SGS_hydro_run_pars !*********************************************************************** subroutine write_SGS_hydro_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=SGS_hydro_run_pars) ! endsubroutine write_SGS_hydro_run_pars !*********************************************************************** subroutine rprint_SGS_hydro(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_nu_LES=0; idiag_Sij2m=0 idiag_fSGSm=0 idiag_fSGSrmsx=0 endif ! ! iname runs through all possible names that may be listed in print.in ! if (lroot.and.ip<1400) print*,'rprint_SGS_hydro: run through parse list' do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'fSGSm',idiag_fSGSm) !call parse_name(iname,cname(iname),cform(iname),'fSGSrmsx',idiag_fSGSrmsx) enddo ! call keep_compiler_quiet(lreset,lwrite) ! endsubroutine rprint_SGS_hydro !*********************************************************************** subroutine pencil_criteria_SGS_hydro ! ! All pencils that the Viscosity module depends on are specified here. ! ! 20-11-04/anders: coded ! lpenc_requested(i_rho)=.true. endsubroutine pencil_criteria_SGS_hydro !*********************************************************************** subroutine pencil_interdep_SGS_hydro(lpencil_in) ! ! Interdependency among pencils from the Viscosity module is specified here. ! ! 20-11-04/anders: coded ! logical, dimension (npencils) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_SGS_hydro !*********************************************************************** subroutine calc_SGS_hydro_force(f,df,p) ! ! Calculate SGS force and store in pencil. ! Most basic pencils should come first, as others may depend on them. ! ! 20-11-04/anders: coded ! use Sub, only: div, div_other, dot2, traceless_strain use Diagnostics, only: sum_mn_name real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! intent(inout) :: f,p intent(out) :: df real, dimension(nx,3) :: tmp1 integer :: j ! ! Divergence of tensor tau. ! Correct only for Cartesian coordinates. ! call div(f,itauSGSRey,p%SGS_force(:,1)) call div_other(f(:,:,:,(/itauSGSRey+1,itauSGSRey+3,itauSGSRey+4/)),p%SGS_force(:,2)) call div_other(f(:,:,:,(/itauSGSRey+2,itauSGSRey+4,itauSGSRey+5/)),p%SGS_force(:,3)) if (lmagnetic) then call div(f,itauSGSMax,tmp1(:,1)) call div_other(f(:,:,:,(/itauSGSMax+1,itauSGSMax+3,itauSGSMax+4/)),tmp1(:,2)) call div_other(f(:,:,:,(/itauSGSMax+2,itauSGSMax+4,itauSGSMax+5/)),tmp1(:,3)) p%SGS_force=p%SGS_force+tmp1 endif if (lSGS_forc_as_aux) f(l1:l2,m,n,iSGS_force:iSGS_force+2)=-p%SGS_force df(l1:l2,m,n,iux:iuz) = df(l1:l2,m,n,iux:iuz) - p%SGS_force ! ! define SGS_heat ! if (lpencil(i_SGS_heat)) p%SGS_heat=0.0 ! call calc_diagnostics_SGS_hydro(p) ! endsubroutine calc_SGS_hydro_force !*********************************************************************** subroutine calc_diagnostics_SGS_hydro(p) use Diagnostics use Sub, only: dot2 type (pencil_case) :: p real, dimension(nx) :: dtp if (idiag_fSGSm/=0) then call dot2(p%SGS_force,dtp) call sum_mn_name(dtp,idiag_fSGSm,lsqrt=.true.) endif endsubroutine calc_diagnostics_SGS_hydro !*********************************************************************** subroutine SGS_hydro_after_boundary(f) use Boundcond, only: update_ghosts, boundconds_x, boundconds_y, boundconds_z use Mpicomm, only: finalize_isendrcv_bdry, initiate_isendrcv_bdry use Sub, only: traceless_strain, div, gij, gij_etc, smooth_kernel real, dimension (mx,my,mz,mfarray) :: f integer :: imn,j real, dimension(nx) :: mij2, penc, sij2 real, dimension(nx,3) :: uu real, dimension(nx,3,3) :: mij, uij, sij logical :: lcommunicate ! ! Requires communication to be finished! ! !print*, 'in SGS_hydro_after_boundary' ! call boundconds_x(f,iux,iuz) ! call initiate_isendrcv_bdry(f,iux,iuz) ! lcommunicate=.true. do imn=1,nyz ! n = nn(imn) m = mm(imn) ! if (lcommunicate) then if (necessary(imn)) then call finalize_isendrcv_bdry(f,iux,iuz) call boundconds_y(f,iuy,iuy) call boundconds_z(f,iuz,iuz) lcommunicate=.false. endif endif do j=iux,iuz call smooth_kernel(f,j,f(l1:l2,m,n,iuusmooth+j-iux),smth_kernel) enddo if (lmagnetic) then do j=iax,iaz call smooth_kernel(f,j,f(l1:l2,m,n,iaasmooth+j-iax),smth_kernel) enddo endif enddo if (lmagnetic) then if (iaasmooth==iuusmooth+3) then call update_ghosts(f,iuusmooth,iaasmooth+2) else call update_ghosts(f,iuusmooth,iuusmooth+2) call update_ghosts(f,iaasmooth,iaasmooth+2) endif else call update_ghosts(f,iuusmooth,iuusmooth+2) endif do n=n1,n2; do m=m1,m2 call gij(f,iuusmooth,uij,1) call div(f,iuusmooth,penc) uu=f(l1:l2,m,n,iuusmooth:iuusmooth+2) call traceless_strain(uij,penc,sij,uu)!,lshear_rateofstrain) sij2=sum(sum(sij**2,2),2) call tauij_SGS(uij,sij2,cReyStress,f,itauSGSRey) ! remember rho if (lmagnetic) then call gij_etc(f,iaasmooth,BIJ=uij) call traceless_strain(uij,sij=mij) mij2=sum(sum(mij**2,2),2) call tauij_SGS(uij,mij2,cMaxStress,f,itauSGSMax) endif enddo; enddo if (lmagnetic) then if (itauSGSMax==itauSGSRey+6) then call update_ghosts(f,itauSGSRey,itauSGSMax+5) else call update_ghosts(f,itauSGSRey,itauSGSRey+5) call update_ghosts(f,itauSGSMax,itauSGSMax+5) endif else call update_ghosts(f,itauSGSRey,itauSGSRey+5) endif endsubroutine SGS_hydro_after_boundary !*********************************************************************** subroutine tauij_SGS(uij,sij2,coef,f,k,rho) real, dimension(nx,3,3) :: uij real, dimension(nx) :: sij2 real, dimension(nx), optional :: rho real, dimension (mx,my,mz,mfarray) :: f real :: coef integer :: k real, dimension(nx) :: E_SGS,uij2 E_SGS=coef*max(dx,dy,dz)**2*sqrt(2.*sij2) if (present(rho)) E_SGS=E_SGS*rho uij2=sum(sum(uij**2,2),2) where (uij2==0) uij2=1. f(l1:l2,m,n,k ) = E_SGS*(sum(uij(:,1,:)**2,2) /uij2 - 1./3.) ! tau(1,1) f(l1:l2,m,n,k+1) = E_SGS*(sum(uij(:,1,:)*uij(:,2,:),2)/uij2 ) ! tau(1,2) f(l1:l2,m,n,k+2) = E_SGS*(sum(uij(:,1,:)*uij(:,3,:),2)/uij2 ) ! tau(1,3) f(l1:l2,m,n,k+3) = E_SGS*(sum(uij(:,2,:)**2,2) /uij2 - 1./3.) ! tau(2,2) f(l1:l2,m,n,k+4) = E_SGS*(sum(uij(:,2,:)*uij(:,3,:),2)/uij2 ) ! tau(2,3) f(l1:l2,m,n,k+5) = E_SGS*(sum(uij(:,3,:)**2,2) /uij2 - 1./3.) ! tau(3,3) endsubroutine tauij_SGS !*********************************************************************** subroutine SGS_hydro_before_boundary(f) real, dimension (mx,my,mz,mfarray) :: f call keep_compiler_quiet(f) ! endsubroutine SGS_hydro_before_boundary !*********************************************************************** subroutine calc_SGS_hydro_heat(df,p,Hmax) ! real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p ! real, dimension (nx) :: Hmax ! intent(in) :: df,p,Hmax ! call keep_compiler_quiet(df) call keep_compiler_quiet(p) call keep_compiler_quiet(Hmax) ! endsubroutine calc_SGS_hydro_heat !*********************************************************************** endmodule SGS_hydro