! $Id$ ! ! This module applies a sixth order hyperviscosity to the equation ! of motion (following Haugen & Brandenburg 2004). This hyperviscosity ! ensures that the energy dissipation rate is positive define everywhere. ! ! The rate of strain tensor ! S^(3) = (-nab^2)^2*S ! is a high order generalisation of the first order operator ! 2*S_ij = u_i,j + u_j,i - 2/3*delta_ij*div(u) ! ! Derivatives are taken in Fourier space. ! !** 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 :: lhyperviscosity_strict=.true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 3 ! !*************************************************************** module Hypervisc_strict ! use Cdata use Cparam use Fourier use Messages ! implicit none ! include 'hypervisc_strict.h' ! contains !*********************************************************************** subroutine register_hypervisc_strict() ! ! Set up indices for hyperviscosity auxiliary slots. ! ! 20-aug-07/anders: coded ! use FArrayManager ! if (lroot) call svn_id( & "$Id$") ! ! Set indices for auxiliary variables. ! call farray_register_auxiliary('hypvis',ihypvis,vector=3) ! endsubroutine register_hypervisc_strict !*********************************************************************** subroutine hyperviscosity_strict(f) ! ! Apply momentum-conserving, symmetric, sixth order hyperviscosity with ! positive define heating rate (see Haugen & Brandenburg 2004). ! ! 20-aug-2007/anders: coded ! use Fourier use Sub ! real, dimension (mx,my,mz,mfarray) :: f ! real, dimension (nx,ny,nz) :: uxhat_re, uxhat_im real, dimension (nx,ny,nz) :: uyhat_re, uyhat_im real, dimension (nx,ny,nz) :: uzhat_re, uzhat_im real, dimension (nx,ny,nz) :: duxhat_re, duxhat_im real, dimension (nx,ny,nz) :: duyhat_re, duyhat_im real, dimension (nx,ny,nz) :: duzhat_re, duzhat_im real, dimension (nx,ny,nz) :: tmp real :: kx, ky, kz complex :: del2del2del2, del2del2divu integer :: ikx, iky, ikz real, dimension (nx,3) :: del6u ! ! Identify version ! if (lroot .and. ip<10) call svn_id( & "$Id$") ! ! Derivatives are taken in k-space due to the complicated cross terms. ! uxhat_re=f(l1:l2,m1:m2,n1:n2,iuu ) uxhat_im=0.0 uyhat_re=f(l1:l2,m1:m2,n1:n2,iuu+1) uyhat_im=0.0 uzhat_re=f(l1:l2,m1:m2,n1:n2,iuu+2) uzhat_im=0.0 ! if (lshear) then call fourier_transform_shear(uxhat_re,uxhat_im) call fourier_transform_shear(uyhat_re,uyhat_im) call fourier_transform_shear(uzhat_re,uzhat_im) else call fourier_transform(uxhat_re,uxhat_im) call fourier_transform(uyhat_re,uyhat_im) call fourier_transform(uzhat_re,uzhat_im) endif ! ! Construct hyperviscous acceleration ! ! f_visc = mu3/rho*(del2(del2(del2(u)))+1/3*del2(del2(grad(div(u)))) ! do ikz=1,nz; do iky=1,ny; do ikx=1,nx ! if (lshear) then kx=kx_fft(ikx)+deltay/Lx*ky_fft(iky+ipy*ny) ky=ky_fft(iky+ipy*ny) kz=kz_fft(ikz) else kx=kx_fft(ikx) ky=ky_fft(iky+ipy*ny) kz=kz_fft(ikz+ipz*nz) endif ! if (lshear) then ! if (nzgrid/=1) then ! 3-D with shear, FFT order (kz,ky',kx) ! ! del2(del2(del2(u))) ! del2del2del2=-(kx**2+ky**2+kz**2)**3 ! duxhat_re(ikz,iky,ikx)=del2del2del2*uxhat_re(ikz,iky,ikx) duxhat_im(ikz,iky,ikx)=del2del2del2*uxhat_im(ikz,iky,ikx) duyhat_re(ikz,iky,ikx)=del2del2del2*uyhat_re(ikz,iky,ikx) duyhat_im(ikz,iky,ikx)=del2del2del2*uyhat_im(ikz,iky,ikx) duzhat_re(ikz,iky,ikx)=del2del2del2*uzhat_re(ikz,iky,ikx) duzhat_im(ikz,iky,ikx)=del2del2del2*uzhat_im(ikz,iky,ikx) ! ! del2(del2(grad(div(u)))) ! del2del2divu= -1/3.* & ! i*1/3*del2(del2(divu)) operator (kx**2+ky**2+kz**2)**2* & (kx*cmplx(uxhat_re(ikz,iky,ikx),uxhat_im(ikz,iky,ikx)) + & ky*cmplx(uyhat_re(ikz,iky,ikx),uyhat_im(ikz,iky,ikx)) + & kz*cmplx(uzhat_re(ikz,iky,ikx),uzhat_im(ikz,iky,ikx)) ) ! duxhat_re(ikz,iky,ikx)=duxhat_re(ikz,iky,ikx)+ real(kx*del2del2divu) duxhat_im(ikz,iky,ikx)=duxhat_im(ikz,iky,ikx)+aimag(kx*del2del2divu) duyhat_re(ikz,iky,ikx)=duyhat_re(ikz,iky,ikx)+ real(ky*del2del2divu) duyhat_im(ikz,iky,ikx)=duyhat_im(ikz,iky,ikx)+aimag(ky*del2del2divu) duzhat_re(ikz,iky,ikx)=duzhat_re(ikz,iky,ikx)+ real(kz*del2del2divu) duzhat_im(ikz,iky,ikx)=duzhat_im(ikz,iky,ikx)+aimag(kz*del2del2divu) ! else ! 2-D with shear, FFT order (kx,ky',kz) ! ! del2(del2(del2(u))) ! del2del2del2=-(kx**2+ky**2+kz**2)**3 ! duxhat_re(ikx,iky,ikz)=del2del2del2*uxhat_re(ikx,iky,ikz) duxhat_im(ikx,iky,ikz)=del2del2del2*uxhat_im(ikx,iky,ikz) duyhat_re(ikx,iky,ikz)=del2del2del2*uyhat_re(ikx,iky,ikz) duyhat_im(ikx,iky,ikz)=del2del2del2*uyhat_im(ikx,iky,ikz) duzhat_re(ikx,iky,ikz)=del2del2del2*uzhat_re(ikx,iky,ikz) duzhat_im(ikx,iky,ikz)=del2del2del2*uzhat_im(ikx,iky,ikz) ! ! del2(del2(grad(div(u)))) ! del2del2divu= -1/3.* & ! i*1/3*del2(del2(divu)) operator (kx**2+ky**2+kz**2)**2* & (kx*cmplx(uxhat_re(ikx,iky,ikz),uxhat_im(ikx,iky,ikz)) + & ky*cmplx(uyhat_re(ikx,iky,ikz),uyhat_im(ikx,iky,ikz)) + & kz*cmplx(uzhat_re(ikx,iky,ikz),uzhat_im(ikx,iky,ikz)) ) ! duxhat_re(ikx,iky,ikz)=duxhat_re(ikx,iky,ikz)+ real(kx*del2del2divu) duxhat_im(ikx,iky,ikz)=duxhat_im(ikx,iky,ikz)+aimag(kx*del2del2divu) duyhat_re(ikx,iky,ikz)=duyhat_re(ikx,iky,ikz)+ real(ky*del2del2divu) duyhat_im(ikx,iky,ikz)=duyhat_im(ikx,iky,ikz)+aimag(ky*del2del2divu) duzhat_re(ikx,iky,ikz)=duzhat_re(ikx,iky,ikz)+ real(kz*del2del2divu) duzhat_im(ikx,iky,ikz)=duzhat_im(ikx,iky,ikz)+aimag(kz*del2del2divu) ! endif ! else ! No shear ! if (nzgrid/=1) then !3D - FFT has put in (kz,kx,ky) order. ! ! del2(del2(del2(u))) ! del2del2del2=-(kx**2+ky**2+kz**2)**3 ! duxhat_re(ikz,ikx,iky)=del2del2del2*uxhat_re(ikz,ikx,iky) duxhat_im(ikz,ikx,iky)=del2del2del2*uxhat_im(ikz,ikx,iky) duyhat_re(ikz,ikx,iky)=del2del2del2*uyhat_re(ikz,ikx,iky) duyhat_im(ikz,ikx,iky)=del2del2del2*uyhat_im(ikz,ikx,iky) duzhat_re(ikz,ikx,iky)=del2del2del2*uzhat_re(ikz,ikx,iky) duzhat_im(ikz,ikx,iky)=del2del2del2*uzhat_im(ikz,ikx,iky) ! ! del2(del2(grad(div(u)))) ! del2del2divu= -1/3.* & ! i*1/3*del2(del2(divu)) operator (kx**2+ky**2+kz**2)**2* & (kx*cmplx(uxhat_re(ikz,ikx,iky),uxhat_im(ikz,ikx,iky)) + & ky*cmplx(uyhat_re(ikz,ikx,iky),uyhat_im(ikz,ikx,iky)) + & kz*cmplx(uzhat_re(ikz,ikx,iky),uzhat_im(ikz,ikx,iky)) ) ! duxhat_re(ikz,ikx,iky)=duxhat_re(ikz,ikx,iky)+ real(kx*del2del2divu) duxhat_im(ikz,ikx,iky)=duxhat_im(ikz,ikx,iky)+aimag(kx*del2del2divu) duyhat_re(ikz,ikx,iky)=duyhat_re(ikz,ikx,iky)+ real(ky*del2del2divu) duyhat_im(ikz,ikx,iky)=duyhat_im(ikz,ikx,iky)+aimag(ky*del2del2divu) duzhat_re(ikz,ikx,iky)=duzhat_re(ikz,ikx,iky)+ real(kz*del2del2divu) duzhat_im(ikz,ikx,iky)=duzhat_im(ikz,ikx,iky)+aimag(kz*del2del2divu) ! else !2D - order (kx,ky,kz) order ! ! del2(del2(del2(u))) ! del2del2del2=-(kx**2+ky**2+kz**2)**3 ! duxhat_re(ikx,iky,ikz)=del2del2del2*uxhat_re(ikx,iky,ikz) duxhat_im(ikx,iky,ikz)=del2del2del2*uxhat_im(ikx,iky,ikz) duyhat_re(ikx,iky,ikz)=del2del2del2*uyhat_re(ikx,iky,ikz) duyhat_im(ikx,iky,ikz)=del2del2del2*uyhat_im(ikx,iky,ikz) duzhat_re(ikx,iky,ikz)=del2del2del2*uzhat_re(ikx,iky,ikz) duzhat_im(ikx,iky,ikz)=del2del2del2*uzhat_im(ikx,iky,ikz) ! ! del2(del2(grad(div(u)))) ! del2del2divu= -1/3.* & ! i*1/3*del2(del2(divu)) operator (kx**2+ky**2+kz**2)**2* & (kx*cmplx(uxhat_re(ikx,iky,ikz),uxhat_im(ikx,iky,ikz)) + & ky*cmplx(uyhat_re(ikx,iky,ikz),uyhat_im(ikx,iky,ikz)) + & kz*cmplx(uzhat_re(ikx,iky,ikz),uzhat_im(ikx,iky,ikz)) ) ! duxhat_re(ikx,iky,ikz)=duxhat_re(ikx,iky,ikz)+ real(kx*del2del2divu) duxhat_im(ikx,iky,ikz)=duxhat_im(ikx,iky,ikz)+aimag(kx*del2del2divu) duyhat_re(ikx,iky,ikz)=duyhat_re(ikx,iky,ikz)+ real(ky*del2del2divu) duyhat_im(ikx,iky,ikz)=duyhat_im(ikx,iky,ikz)+aimag(ky*del2del2divu) duzhat_re(ikx,iky,ikz)=duzhat_re(ikx,iky,ikz)+ real(kz*del2del2divu) duzhat_im(ikx,iky,ikz)=duzhat_im(ikx,iky,ikz)+aimag(kz*del2del2divu) ! endif ! endif ! enddo; enddo; enddo ! ! Inverse transform (to real space). The Viscosity module multiplies the ! result with mu3/rho. ! if (lshear) then call fourier_transform_shear(duxhat_re,duxhat_im,linv=.true.) f(l1:l2,m1:m2,n1:n2,ihypvis)=duxhat_re call fourier_transform_shear(duyhat_re,duyhat_im,linv=.true.) f(l1:l2,m1:m2,n1:n2,ihypvis+1)=duyhat_re call fourier_transform_shear(duzhat_re,duzhat_im,linv=.true.) f(l1:l2,m1:m2,n1:n2,ihypvis+2)=duzhat_re else call fourier_transform(duxhat_re,duxhat_im,linv=.true.) f(l1:l2,m1:m2,n1:n2,ihypvis)=duxhat_re call fourier_transform(duyhat_re,duyhat_im,linv=.true.) f(l1:l2,m1:m2,n1:n2,ihypvis+1)=duyhat_re call fourier_transform(duzhat_re,duzhat_im,linv=.true.) f(l1:l2,m1:m2,n1:n2,ihypvis+2)=duzhat_re endif ! ! The imaginary part of the functions should be zero once transformed ! back to real space, since the differential operators preserve ! complex conjugarity of the scales (kx,ky,kz) and (-kx,-ky,-kz) ! if (ip<=10) then print'(A,4e15.5)', & 'hypervisc_strict: min(duxhat_re), max(duxhat_re), '// & 'min(duxhat_im), max(duxhat_im)', & minval(duxhat_re), maxval(duxhat_re), & minval(duxhat_im), maxval(duxhat_im) print'(A,4e15.5)', & 'hypervisc_strict: min(duyhat_re), max(duyhat_re), '// & 'min(duyhat_im), max(duyhat_im)', & minval(duyhat_re), maxval(duyhat_re), & minval(duyhat_im), maxval(duyhat_im) print'(A,4e15.5)', & 'hypervisc_strict: min(duzhat_re), max(duzhat_re), '// & 'min(duzhat_im), max(duzhat_im)', & minval(duzhat_re), maxval(duzhat_re), & minval(duzhat_im), maxval(duzhat_im) endif ! endsubroutine hyperviscosity_strict !*********************************************************************** endmodule Hypervisc_strict