mod_viscosity.t

Go to the documentation of this file.

!>   The module add viscous source terms and check time step
!>
!>   Viscous forces in the momentum equations:
!>   d m_i/dt +=  - div (vc_mu * PI)
!>   !! Viscous work in the energy equation:
!>   !! de/dt    += - div (v . vc_mu * PI)
!>   where the PI stress tensor is
!>   PI_i,j = - (dv_j/dx_i + dv_i/dx_j) + (2/3)*Sum_k dv_k/dx_k
!>   where vc_mu is the dynamic viscosity coefficient (g cm^-1 s^-1).
module mod_viscosity
  use mod_comm_lib, only: mpistop
  implicit none
 
  !> Viscosity coefficient
  double precision, public :: vc_mu = 1.d0
 
  !> Index of the density (in the w array)
  integer, private, parameter              :: rho_ = 1
 
  !> Indices of the momentum density
  integer, allocatable, private, protected :: mom(:)
 
  !> Index of the energy density (-1 if not present)
  integer, private, protected              :: e_
 
  !> fourth order
  logical :: vc_4th_order = .false.
 
  !> source split or not
  logical :: vc_split= .false.
 
  !> whether to compute the viscous terms as
  !> fluxes (ie in the div on the LHS), or not (by default)
  logical :: viscindiv= .false.
 
 
  ! Public methods
  public :: visc_get_flux_prim
 
contains
  !> Read this module"s parameters from a file
  subroutine vc_params_read(files)
    use mod_global_parameters, only: unitpar
    character(len=*), intent(in) :: files(:)
    integer                      :: n
 
    namelist /vc_list/ vc_mu, vc_4th_order, vc_split, viscindiv
 
    do n = 1, size(files)
       open(unitpar, file=trim(files(n)), status="old")
       read(unitpar, vc_list, end=111)
111    close(unitpar)
    end do
 
  end subroutine vc_params_read
 
  !> Initialize the module
  subroutine viscosity_init(phys_wider_stencil)
    use mod_global_parameters
    integer, intent(inout) :: phys_wider_stencil
    integer :: nwx,idir
 
    call vc_params_read(par_files)
 
    if(vc_split) any_source_split=.true.
 
    ! Determine flux variables
    nwx = 1                  ! rho (density)
 
    allocate(mom(ndir))
    do idir = 1, ndir
       nwx    = nwx + 1
       mom(idir) = nwx       ! momentum density
    end do
 
    nwx = nwx + 1
    e_     = nwx          ! energy density
 
    if (viscindiv) then
      ! to compute the derivatives from left and right upwinded values
      phys_wider_stencil = 1
    end if
 
  end subroutine viscosity_init
 
  subroutine viscosity_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
       energy,qsourcesplit,active)
  ! Add viscosity source in isotropic Newtonian fluids to w within ixO
  ! neglecting bulk viscosity
  ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
    use mod_global_parameters
    use mod_geometry
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
    double precision, intent(inout) :: w(ixI^S,1:nw)
    logical, intent(in) :: energy,qsourcesplit
    logical, intent(inout) :: active
 
    double precision:: lambda(ixI^S,ndir,ndir),tmp(ixI^S),tmp2(ixI^S),v(ixI^S,ndir),vlambda(ixI^S,ndir)
    integer:: ix^L,idim,idir,jdir,iw
 
    if (viscindiv) return
 
    if(qsourcesplit .eqv. vc_split) then
      active = .true.
      ! standard case, textbook viscosity
      ! Calculating viscosity sources
      if(.not.vc_4th_order) then
        ! involves second derivatives, two extra layers
        ix^l=ixo^l^ladd2;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition, 2 layers needed")
        ix^l=ixo^l^ladd1;
      else
        ! involves second derivatives, four extra layers
        ix^l=ixo^l^ladd4;
        if({ iximin^d>ixmin^d .or. iximax^d<ixmax^d|.or.})&
          call mpistop("error for viscous source addition"//&
          "requested fourth order gradients: 4 layers needed")
        ix^l=ixo^l^ladd2;
      end if
 
      ! get velocity
      do idir=1,ndir
        v(ixi^s,idir)=wct(ixi^s,mom(idir))/wct(ixi^s,rho_)
      end do
 
      ! construct lambda tensor: lambda_ij = gradv_ij + gradv_ji
      ! initialize
      lambda=zero
 
      !next construct
      do idim=1,ndim; do idir=1,ndir
      ! Calculate velocity gradient tensor within ixL: gradv= grad v,
      ! thus gradv_ij=d_j v_i
        tmp(ixi^s)=v(ixi^s,idir)
        ! Correction for Christoffel terms in non-cartesian
        if (coordinate==cylindrical .and. idim==r_  .and. idir==phi_  ) tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
        if (coordinate==spherical) then
          if     (idim==r_  .and. (idir==2 .or. idir==phi_)) then
            tmp(ixi^s) = tmp(ixi^s)/x(ixi^s,1)
{^nooned
          elseif (idim==2  .and. idir==phi_) then
            tmp(ixi^s)=tmp(ixi^s)/dsin(x(ixi^s,2))
}
          endif
        endif
        call gradient(tmp,ixi^l,ix^l,idim,tmp2)
        ! Correction for Christoffel terms in non-cartesian
        if (coordinate==cylindrical .and. idim==r_  .and. idir==phi_  ) tmp2(ix^s)=tmp2(ix^s)*x(ix^s,1)
        if (coordinate==cylindrical .and. idim==phi_ .and. idir==phi_ ) tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
        if (coordinate==spherical) then
          if (idim==r_  .and. (idir==2 .or. idir==phi_)) then
            tmp2(ix^s) = tmp2(ix^s)*x(ix^s,1)
{^nooned
          elseif (idim==2  .and. idir==phi_ ) then
            tmp2(ix^s)=tmp2(ix^s)*dsin(x(ix^s,2))
          elseif (idim==2   .and. idir==2   ) then
            tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)
          elseif (idim==phi_.and. idir==phi_) then
            tmp2(ix^s)=tmp2(ix^s)+v(ix^s,r_)/x(ix^s,1)+v(ix^s,2)/(x(ix^s,1)*dtan(x(ix^s,2)))
}
          endif
        endif
        lambda(ix^s,idim,idir)= lambda(ix^s,idim,idir)+ tmp2(ix^s)
        lambda(ix^s,idir,idim)= lambda(ix^s,idir,idim)+ tmp2(ix^s)
      enddo; enddo;
 
      ! Multiply lambda with viscosity coefficient and dt
      lambda(ix^s,1:ndir,1:ndir)=lambda(ix^s,1:ndir,1:ndir)*vc_mu*qdt
 
      !calculate div v term through trace action separately
      ! rq : it is safe to use the trace rather than compute the divergence
      !      since we always retrieve the divergence (even with the
      !      Christoffel terms)
      tmp=0.d0
      do idir=1,ndir
         tmp(ix^s)=tmp(ix^s)+lambda(ix^s,idir,idir)
      end do
      tmp(ix^s)=tmp(ix^s)/3.d0
 
      !substract trace from diagonal elements
      do idir=1,ndir
         lambda(ix^s,idir,idir)=lambda(ix^s,idir,idir)-tmp(ix^s)
      enddo
 
      ! dm/dt= +div(mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v * kr)
      ! hence m_j=m_j+d_i tensor_ji
      do idir=1,ndir
        do idim=1,ndim
              tmp(ix^s)=lambda(ix^s,idir,idim)
              ! Correction for divergence of a tensor
              if (coordinate==cylindrical .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)
              if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_              ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
              if (coordinate==spherical) then
                if (idim==r_ .and. idir==r_                 ) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**two
                if (idim==r_ .and. (idir==2 .or. idir==phi_)) tmp(ix^s) = tmp(ix^s)*x(ix^s,1)**3.d0
{^nooned
                if (idim==2  .and. (idir==r_ .or. idir==2))   tmp(ix^s) = tmp(ix^s)*dsin(x(ix^s,2))
                if (idim==2  .and. idir==phi_               ) tmp(ix^s) = tmp(ix^s)*dsin(x(ix^s,2))**two
}
              endif
              call gradient(tmp,ixi^l,ixo^l,idim,tmp2)
              ! Correction for divergence of a tensor
              if (coordinate==cylindrical .and. idim==r_ .and. (idir==r_ .or. idir==z_)) tmp2(ixo^s) = tmp2(ixo^s)/x(ixo^s,1)
              if (coordinate==cylindrical .and. idim==r_ .and. idir==phi_              ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
              if (coordinate==spherical) then
                if (idim==r_ .and. idir==r_                 ) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**two)
                if (idim==r_ .and. (idir==2 .or. idir==phi_)) tmp2(ixo^s) = tmp2(ixo^s)/(x(ixo^s,1)**3.d0)
{^nooned
                if (idim==2  .and. (idir==r_ .or. idir==2))   tmp2(ixo^s) = tmp2(ixo^s)/(dsin(x(ixo^s,2)))
                if (idim==2  .and. idir==phi_               ) tmp2(ixo^s) = tmp2(ixo^s)/(dsin(x(ixo^s,2))**two)
}
              endif
              w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp2(ixo^s)
        enddo
        ! Correction for geometrical terms in the div of a tensor
        if (coordinate==cylindrical .and. idir==r_  ) w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-lambda(ixo^s,phi_,phi_)/x(ixo^s,1)
        if (coordinate==spherical   .and. idir==r_  ) w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-(lambda(ixo^s,2,2)+lambda(ixo^s,phi_,phi_))/x(ixo^s,1)
{^nooned
        if (coordinate==spherical   .and. idir==2   ) w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-lambda(ixo^s,phi_,phi_)/(x(ixo^s,1)/dtan(x(ixo^s,2)))
}
      end do
 
      if(energy) then
        ! de/dt= +div(v.dot.[mu*[d_j v_i+d_i v_j]-(2*mu/3)* div v *kr])
        ! thus e=e+d_i v_j tensor_ji
        vlambda=0.d0
        do idim=1,ndim
          do idir=1,ndir
             vlambda(ixi^s,idim)=vlambda(ixi^s,idim)+v(ixi^s,idir)*lambda(ixi^s,idir,idim)
          end do
        end do
        call divvector(vlambda,ixi^l,ixo^l,tmp2)
        w(ixo^s,e_)=w(ixo^s,e_)+tmp2(ixo^s)
      end if
    end if
 
  end subroutine viscosity_add_source
 
  subroutine viscosity_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
    ! Check diffusion time limit for dt < dtdiffpar * dx**2 / (mu/rho)
    use mod_global_parameters
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: dx^D, x(ixI^S,1:ndim)
    double precision, intent(in) :: w(ixI^S,1:nw)
    double precision, intent(inout) :: dtnew
 
    double precision :: tmp(ixI^S)
    double precision:: dtdiff_visc, dxinv2(1:ndim)
    integer:: idim
 
    ! Calculate the kinematic viscosity tmp=mu/rho
 
    tmp(ixo^s)=vc_mu/w(ixo^s,rho_)
 
    ^d&dxinv2(^d)=one/dx^d**2;
    do idim=1,ndim
       dtdiff_visc=dtdiffpar/maxval(tmp(ixo^s)*dxinv2(idim))
       ! limit the time step
       dtnew=min(dtnew,dtdiff_visc)
    enddo
 
  end subroutine viscosity_get_dt
 
  ! viscInDiv
  ! Get the viscous stress tensor terms in the idim direction
  ! Beware : a priori, won't work for ndir /= ndim
  ! Rq : we work with primitive w variables here
  ! Rq : ixO^L is already extended by 1 unit in the direction we work on
 
  subroutine visc_get_flux_prim(w, x, ixI^L, ixO^L, idim, f, energy)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixi^l, ixo^l, idim
    double precision, intent(in)    :: w(ixi^s, 1:nw), x(ixi^s, 1:^nd)
    double precision, intent(inout) :: f(ixi^s, nwflux)
    logical, intent(in) :: energy
    integer                         :: idir, i
    double precision :: v(ixi^s,1:ndir)
 
    double precision                :: divergence(ixi^s)
 
    double precision:: lambda(ixi^s,ndir) !, tmp(ixI^S) !gradV(ixI^S,ndir,ndir)
 
    if (.not. viscindiv) return
 
    do i=1,ndir
     v(ixi^s,i)=w(ixi^s,i+1)
    enddo
    call divvector(v,ixi^l,ixo^l,divergence)
 
    call get_crossgrad(ixi^l,ixo^l,x,w,idim,lambda)
    lambda(ixo^s,idim) = lambda(ixo^s,idim) - (2.d0/3.d0) * divergence(ixo^s)
 
    ! Compute the idim-th row of the viscous stress tensor
    do idir = 1, ndir
      f(ixo^s, mom(idir)) = f(ixo^s, mom(idir)) - vc_mu * lambda(ixo^s,idir)
      if (energy) f(ixo^s, e_) = f(ixo^s, e_) - vc_mu * lambda(ixo^s,idir) * v(ixi^s,idir)
    enddo
 
  end subroutine visc_get_flux_prim
 
  ! Compute the cross term ( d_i v_j + d_j v_i in Cartesian BUT NOT IN
  ! CYLINDRICAL AND SPHERICAL )
  subroutine get_crossgrad(ixI^L,ixO^L,x,w,idim,cross)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: w(ixI^S, 1:nw), x(ixI^S, 1:ndim)
    double precision, intent(out)   :: cross(ixI^S,ndir)
 
    double precision :: tmp(ixI^S), v(ixI^S)
    integer :: idir
 
    if (ndir/=ndim) call mpistop("This formula are probably wrong for ndim/=ndir")
    ! Beware also, we work w/ the angle as the 3rd component in cylindrical
    ! and the colatitude as the 2nd one in spherical
    cross(ixi^s,:)=zero
    tmp(ixi^s)=zero
    select case(coordinate)
    case (cartesian,cartesian_stretched)
      call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
    case (cylindrical)
      if (idim==1) then
        ! for rr and rz
        call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
        ! then we overwrite rth w/ the correct expression
        {^nooned
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
        }
      elseif (idim==2) then
        ! thr (idem as above)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        {^nooned
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,1)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
        ! thth
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=two*tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
        !thz
        v(ixi^s)=w(ixi^s,mom(3)) ! v_z
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        }
        cross(ixi^s,3)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))  ! v_th
        {^ifthreed
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
        }
        cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)
        {^ifthreed
      elseif (idim==3) then
        ! for zz and rz
        call cart_cross_grad(ixi^l,ixo^l,x,w,idim,cross)
        ! then we overwrite zth w/ the correct expression
        !thz
        v(ixi^s)=w(ixi^s,mom(3)) ! v_z
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))  ! v_th
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_z
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)
        }
      endif
    case (spherical)
      if (idim==1) then
        ! rr (normal, simply 2 * dr vr)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=two*tmp(ixi^s)
        !rth
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        {^nooned
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*x(ixi^s,1)
        }
        {^ifthreed
        !rph
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
        cross(ixi^s,3)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*x(ixi^s,1)
        }
      elseif (idim==2) then
        ! thr
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        {^nooned
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
        cross(ixi^s,1)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
        ! thth
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=two*tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        cross(ixi^s,2)=cross(ixi^s,2)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
        }
        {^ifthreed
        !thph
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
        cross(ixi^s,3)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,3)=cross(ixi^s,3)+tmp(ixi^s)*dsin(x(ixi^s,2))
        }
        {^ifthreed
      elseif (idim==3) then
        !phr
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
        cross(ixi^s,1)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
        call gradient(v,ixi^l,ixo^l,1,tmp) ! d_r
        cross(ixi^s,1)=cross(ixi^s,1)+tmp(ixi^s)*x(ixi^s,1)
        !phth
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
        cross(ixi^s,2)=tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
        call gradient(v,ixi^l,ixo^l,2,tmp) ! d_th
        cross(ixi^s,2)=cross(ixi^s,2)+tmp(ixi^s)*dsin(x(ixi^s,2))
        !phph
        v(ixi^s)=w(ixi^s,mom(3)) ! v_ph
        call gradient(v,ixi^l,ixo^l,3,tmp) ! d_phi
        cross(ixi^s,3)=two*tmp(ixi^s)
        v(ixi^s)=w(ixi^s,mom(1)) ! v_r
        cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/x(ixi^s,1) ! + 2 vr/r
        v(ixi^s)=w(ixi^s,mom(2)) ! v_th
        cross(ixi^s,3)=cross(ixi^s,3)+two*v(ixi^s)/(x(ixi^s,1)*dtan(x(ixi^s,2))) ! + 2 vth/(rtan(th))
        }
      endif
    case default
      call mpistop("Unknown geometry specified")
    end select
 
  end subroutine get_crossgrad
 
  !> yields d_i v_j + d_j v_i for a given i, OK in Cartesian and for some
  !> tensor terms in cylindrical (rr & rz) and in spherical (rr)
  subroutine cart_cross_grad(ixI^L,ixO^L,x,w,idim,cross)
    use mod_global_parameters
    use mod_geometry
    integer, intent(in)             :: ixI^L, ixO^L, idim
    double precision, intent(in)    :: w(ixI^S, 1:nw), x(ixI^S, 1:^ND)
    double precision, intent(out)   :: cross(ixI^S,ndir)
 
    double precision :: tmp(ixI^S), v(ixI^S)
    integer :: idir
 
    v(ixi^s)=w(ixi^s,mom(idim))
    do idir=1,ndir
      call gradient(v,ixi^l,ixo^l,idir,tmp)
      cross(ixo^s,idir)=tmp(ixo^s)
    enddo
    do idir=1,ndir
      v(ixi^s)=w(ixi^s,mom(idir))
      call gradient(v,ixi^l,ixo^l,idim,tmp)
      cross(ixo^s,idir)=cross(ixo^s,idir)+tmp(ixo^s)
    enddo
 
  end subroutine cart_cross_grad
 
  subroutine visc_add_source_geom(qdt, ixI^L, ixO^L, wCT, w, x)
    use mod_global_parameters
    use mod_geometry
    ! w and wCT conservative variables here
    integer, intent(in)             :: ixI^L, ixO^L
    double precision, intent(in)    :: qdt, x(ixI^S, 1:ndim)
    double precision, intent(inout) :: wCT(ixI^S, 1:nw), w(ixI^S, 1:nw)
    ! to change and to set as a parameter in the parfile once the possibility to
    ! solve the equations in an angular momentum conserving form has been
    ! implemented (change tvdlf.t eg)
    double precision :: vv(ixI^S), divergence(ixI^S)
    double precision :: tmp(ixI^S),tmp1(ixI^S)
    integer          :: i
 
    if (.not. viscindiv) return
 
    select case (coordinate)
    case (cylindrical)
      ! thth tensor term - - -
        ! 1st the cross grad term
{^nooned
      call gradient(wct(ixi^s,mom(2)),ixi^l,ixo^l,2,tmp1) ! d_th
      tmp(ixo^s)=two*(tmp1(ixo^s)+wct(ixo^s,mom(1))/x(ixo^s,1)) ! 2 ( d_th v_th / r + vr/r )
        ! 2nd the divergence
      call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divergence)
      tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
      ! s[mr]=-thth/radius
      w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! rth tensor term - - -
      call gradient(wct(ixi^s,mom(1)),ixi^l,ixo^l,2,tmp) ! d_th
      vv(ixi^s)=wct(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
      call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
      ! s[mphi]=+rth/radius
      w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
}
    case (spherical)
      ! thth tensor term - - -
      ! 1st the cross grad term
{^nooned
      call gradient(wct(ixi^s,mom(2)),ixi^l,ixo^l,2,tmp1) ! d_th
      tmp(ixo^s)=two*(tmp1(ixo^s)+wct(ixo^s,mom(1))/x(ixo^s,1)) ! 2 ( 1/r * d_th v_th + vr/r )
      ! 2nd the divergence
      call divvector(wct(ixi^s,mom(1:ndir)),ixi^l,ixo^l,divergence)
      tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
      ! s[mr]=-thth/radius
      w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
}
      ! phiphi tensor term - - -
      ! 1st the cross grad term
{^ifthreed
      call gradient(wct(ixi^s,mom(3)),ixi^l,ixo^l,3,tmp1) ! d_phi
      tmp(ixo^s)=two*(tmp1(ixo^s)+wct(ixo^s,mom(1))/x(ixo^s,1)+wct(ixo^s,mom(2))/(x(ixo^s,1)*dtan(x(ixo^s,2)))) ! 2 ( 1/rsinth * d_ph v_ph + vr/r + vth/rtanth )
}
      ! 2nd the divergence
      tmp(ixo^s) = tmp(ixo^s) - (2.d0/3.d0) * divergence(ixo^s)
      ! s[mr]=-phiphi/radius
      w(ixo^s,mom(1))=w(ixo^s,mom(1))-qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! s[mth]=-cotanth*phiphi/radius
{^nooned
      w(ixo^s,mom(2))=w(ixo^s,mom(2))-qdt*vc_mu*tmp(ixo^s)/(x(ixo^s,1)*dtan(x(ixo^s,2)))
}
      ! rth tensor term - - -
      call gradient(wct(ixi^s,mom(1)),ixi^l,ixo^l,2,tmp) ! d_th (rq : already contains 1/r)
      vv(ixi^s)=wct(ixi^s,mom(2))/x(ixi^s,1)  ! v_th / r
      call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
      ! s[mth]=+rth/radius
      w(ixo^s,mom(2))=w(ixo^s,mom(2))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! rphi tensor term - - -
{^ifthreed
      call gradient(wct(ixi^s,mom(1)),ixi^l,ixo^l,3,tmp) ! d_phi (rq : contains 1/rsin(th))
}
      vv(ixi^s)=wct(ixi^s,mom(3))/x(ixi^s,1) ! v_phi / r
      call gradient(vv,ixi^l,ixo^l,1,tmp1) ! d_r
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*x(ixo^s,1)
      ! s[mphi]=+rphi/radius
      w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/x(ixo^s,1)
      ! phith tensor term - - -
{^ifthreed
      call gradient(wct(ixi^s,mom(2)),ixi^l,ixo^l,3,tmp) ! d_phi
}
{^nooned
      vv(ixi^s)=wct(ixi^s,mom(3))/dsin(x(ixi^s,2)) ! v_ph / sin(th)
      call gradient(vv,ixi^l,ixo^l,2,tmp1) ! d_th
      tmp(ixo^s)=tmp(ixo^s)+tmp1(ixo^s)*dsin(x(ixo^s,2))
      ! s[mphi]=+cotanth*phith/radius
      w(ixo^s,mom(3))=w(ixo^s,mom(3))+qdt*vc_mu*tmp(ixo^s)/(x(ixo^s,1)*dtan(x(ixo^s,2)))
}
    end select
 
  end subroutine visc_add_source_geom
 
end module mod_viscosity