sources_viscosity.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: sources_viscosity.f90                                             #
!#                                                                           #
!# Copyright (C) 2008-2019                                                   #
!# Bjoern Sperling <sperling@astrophysik.uni-kiel.de>                        #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!#                                                                           #
!# This program is free software; you can redistribute it and/or modify      #
!# it under the terms of the GNU General Public License as published by      #
!# the Free Software Foundation; either version 2 of the License, or (at     #
!# your option) any later version.                                           #
!#                                                                           #
!# This program is distributed in the hope that it will be useful, but       #
!# WITHOUT ANY WARRANTY; without even the implied warranty of                #
!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or        #
!# NON INFRINGEMENT.  See the GNU General Public License for more            #
!# details.                                                                  #
!#                                                                           #
!# You should have received a copy of the GNU General Public License         #
!# along with Mesh program; if not, write to the Free Software               #
!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################
!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
MODULE sources_viscosity_mod
  USE sources_base_mod
  USE physics_base_mod
  USE fluxes_base_mod
  USE mesh_base_mod
  USE marray_base_mod
  USE marray_compound_mod
  USE logging_base_mod
  USE common_dict
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  PRIVATE
    ! flags for viscosity model
    INTEGER,PARAMETER :: molecular = 1     ! constant viscosity
    INTEGER,PARAMETER :: alpha     = 2     ! Shakura-Sunyaev prescription
    INTEGER,PARAMETER :: beta      = 3     ! Duschl prescription
    INTEGER,PARAMETER :: powerlaw  = 4     ! power law depending on spec. angular momentum
    INTEGER,PARAMETER :: alpha_alt = 5     ! alternative Shakura-Sunyaev
    CHARACTER(LEN=32),PARAMETER, DIMENSION(5) :: viscosity_name = (/ &
                                     "constant viscosity              ", &
                                     "turbulent Shakura-Sunyaev       ", &
                                     "turbulent Duschl                ", &
                                     "power law viscosity             ", &
                                     "alternative Shakura-Sunyaev     "/)


  TYPE, EXTENDS(sources_base) :: sources_viscosity
    CHARACTER(LEN=32) :: source_name = "viscosity of Newtonian fluid"
    CLASS(logging_base), ALLOCATABLE :: viscosity
    CLASS(marray_base), ALLOCATABLE  :: ephir, &    !< azimuthal unit vector */ distance to axis
                                        ephir_tmp, &
                                        dynvis, &   !< dynamic viscosity
                                        kinvis, &   !< kinematic viscosity
                                        bulkvis
    REAL :: dynconst,bulkconst,power
    REAL, DIMENSION(:,:,:), POINTER :: &            !< comp. of stress tensor
          btxx,btyy,btzz,btxy,btxz,btyz,tmp,tmp2,tmp3
    REAL, DIMENSION(:,:), POINTER   :: &
          sxx,syy,szz,sxy,sxz,syz
  CONTAINS
    PROCEDURE :: initsources_viscosity
    PROCEDURE :: infosources
    PROCEDURE :: setoutput
    PROCEDURE :: updateviscosity
    PROCEDURE :: externalsources_single
    PROCEDURE :: calctimestep_single

    PROCEDURE :: finalize
  END TYPE
  !--------------------------------------------------------------------------!
 PUBLIC :: &
       ! types
       sources_viscosity, &
       ! constants
       molecular,alpha,beta,powerlaw,alpha_alt,viscosity_name
  !--------------------------------------------------------------------------!

CONTAINS

  SUBROUTINE initsources_viscosity(this,Mesh,Physics,Fluxes,config,IO)
    USE geometry_base_mod
    USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
    USE geometry_logcylindrical_mod, ONLY: geometry_logcylindrical
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Sources_viscosity) :: this
    CLASS(Mesh_base)    :: Mesh
    CLASS(Physics_base) :: Physics
    CLASS(Fluxes_base)  :: Fluxes
    TYPE(Dict_TYP),POINTER :: config,IO
    INTEGER           :: stype
    !------------------------------------------------------------------------!
    INTEGER           :: err, viscosity_number ,k,l
    !------------------------------------------------------------------------!
    INTENT(IN)        :: mesh,physics,fluxes
    !------------------------------------------------------------------------!
    CALL getattr(config, "stype", stype)
    CALL this%InitLogging(stype,this%source_name)

    SELECT TYPE(phys => physics)
    CLASS IS(physics_eulerisotherm)
      IF (phys%VDIM.EQ.1) &
        CALL this%Error("InitSources_viscosity","viscosity is currently not supported in 1D physics")
    CLASS DEFAULT
      CALL this%Error("InitSources_viscosity", &
        "viscosity is currently only supported in eulerisotherm/euler physics")
    END SELECT

    ! viscosity model
    CALL getattr(config, "vismodel", viscosity_number)
    ALLOCATE(logging_base::this%viscosity)
    CALL this%viscosity%InitLogging(viscosity_number,viscosity_name(viscosity_number))

    ! dynamic viscosity constant
    CALL getattr(config, "dynconst", this%dynconst, 0.1)
    ! bulk viscosity constant
    ! At the moment this only affects the molecular viscosity. The default
    ! value refers to a trace free viscosity tensor.
    ! All other viscosity models should always be trace free. Therefore this
    ! is hard coded and cannot be changed.
    CALL getattr(config, "bulkconst", this%bulkconst, -2./3.*this%dynconst)

    ! power law exponent used for POWERLAW viscosity;
    ! defaults to constant kinematic viscosity, i.e. nu ~ l^0 = const
    CALL getattr(config, "exponent", this%power, 0.0)

    ! set viscous courant number
    CALL getattr(config, "cvis", this%cvis, 0.5)

    ALLOCATE(this%dynvis,this%kinvis,this%bulkvis, &
         this%btxx(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
         this%btyy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
         this%btzz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
         this%btxy(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
         this%btxz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
         this%btyz(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
         stat=err)
    IF (err.NE.0) &
         CALL this%Error("InitSources_viscosity","Memory allocation failed.")

    ! initialize mesh arrays
    this%dynvis  = marray_base()
    this%kinvis  = marray_base()
    this%bulkvis = marray_base()
    this%dynvis%data1d(:)  = this%dynconst
    this%kinvis%data1d(:)  = 0.0
    this%bulkvis%data1d(:) = this%bulkconst

    ! do the model specific initialization
    SELECT CASE(this%viscosity%GetType())
    CASE(molecular)
       ! do nothing
    CASE(alpha,beta,powerlaw)
      ! compute the projector to determine the local angular velocity
      ! ATTENTION: it is always assumed that the fluid rotates around the z-axis
      SELECT TYPE(phys => physics)
      CLASS IS(physics_eulerisotherm)
        ALLOCATE(this%ephir,this%ephir_tmp,stat=err)
        IF (err.NE.0) &
          CALL this%Error("InitSources_viscosity","Memory allocation failed.")
        this%ephir_tmp = marray_base(3)
        this%ephir = marray_base(3)
        ! compute ephi/r (r: distance to z-axis, ephi: azimuthal unit vector)
        ! set cartesian vector components using cartesian coordinates
        this%ephir%data4d(:,:,:,1) = -mesh%cart%bcenter(:,:,:,2) / mesh%radius%bcenter(:,:,:)**2
        this%ephir%data4d(:,:,:,2) = mesh%cart%bcenter(:,:,:,1) / mesh%radius%bcenter(:,:,:)**2
        this%ephir%data4d(:,:,:,3) = 0.0
        ! convert to curvilinear vector components
        CALL mesh%geometry%Convert2Curvilinear(mesh%curv%bcenter,this%ephir%data4d,this%ephir_tmp%data4d)
        ! reduce vector components to match Physics%VDIM
        IF (phys%VDIM.LT.3) THEN
          CALL this%ephir%Destroy()
          this%ephir = marray_base(phys%VDIM)
        END IF
        ! copy vector components actually used
        l = 1
        DO k=1,3
          IF (btest(mesh%VECTOR_COMPONENTS,k-1)) THEN
            this%ephir%data2d(:,l) = this%ephir_tmp%data2d(:,k)
            l = l + 1
          END IF
        END DO
        ! destroy temporary marray
        CALL this%ephir_tmp%Destroy()
        DEALLOCATE(this%ephir_tmp)
      CLASS DEFAULT
        CALL this%Error("InitSources_viscosity",&
               "Physics not supported for alpha-/beta-viscosity")
      END SELECT
    CASE(alpha_alt)
       CALL this%Error("InitSources_viscosity",&
               "alternative alpha-viscosity currently not supported")
!        IF (Physics%VDIM.NE.2) &
!           CALL this%Error("InitSources_viscosity",&
!                "alternative alpha-viscosity works only for flat disks")
   !    IF (.NOT.Timedisc%always_update_bccsound) &
   !       CALL Error(this,"InitSources_viscosity","always_update_bccsound must be enabled in timedisc")
    CASE DEFAULT
       CALL this%Error("InitSources_viscosity",&
            "Viscosity prescription not supported.")
    END SELECT
    ! set initial time to negative value;
    ! this guarantees update of viscosity for time t=0.0
    this%time = -1.0
    ! set stress tensor arrays to zero
    this%btxx(:,:,:) = 0.0
    this%btyy(:,:,:) = 0.0
    this%btzz(:,:,:) = 0.0
    this%btxy(:,:,:) = 0.0
    this%btxz(:,:,:) = 0.0
    this%btyz(:,:,:) = 0.0

    CALL this%SetOutput(mesh,physics,config,io)
    IF (this%GetType().EQ.alpha_alt) this%update_disk_height = .true.

    CALL this%InitSources(mesh,fluxes,physics,config,io)

  END SUBROUTINE initsources_viscosity

  SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Sources_viscosity),INTENT(IN) :: this
    CLASS(Mesh_base) ,INTENT(IN)      :: Mesh
    CLASS(Physics_base),INTENT(IN)    :: Physics
    TYPE(Dict_TYP),POINTER  :: config,IO
    !------------------------------------------------------------------------!
    INTEGER              :: valwrite
    !------------------------------------------------------------------------!
   ! Hier wird der Stresstensor per default in 3D ausgegeben! Vorher nur in 2D plus abhängig
   ! von der Physik auch die Komponenten der dritten Dimension
    CALL getattr(config, "output/stress", valwrite, 0)
    IF (valwrite .EQ. 1) THEN
       CALL setattr(io, "stress:xx", this%btxx(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
       CALL setattr(io, "stress:xy", this%btxy(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
       CALL setattr(io, "stress:xz", this%btxz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
       CALL setattr(io, "stress:yy", this%btyy(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
       CALL setattr(io, "stress:yz", this%btyz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
       CALL setattr(io, "stress:zz", this%btzz(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))
    END IF

    CALL getattr(config, "output/dynvis", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io,"dynvis", this%dynvis%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))

    CALL getattr(config, "output/kinvis", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "kinvis", this%kinvis%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))

    CALL getattr(config, "output/bulkvis", valwrite, 0)
    IF (valwrite .EQ. 1) &
         CALL setattr(io, "bulkvis", this%bulkvis%data3d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX))

  END SUBROUTINE setoutput


  SUBROUTINE infosources(this,Mesh)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Sources_viscosity),INTENT(IN) :: this
    CLASS(Mesh_base),INTENT(IN)         :: Mesh
    !------------------------------------------------------------------------!
    CHARACTER(LEN=32) :: dynconst_str,bulkconst_str
    !------------------------------------------------------------------------!
    WRITE (dynconst_str,'(ES9.2)') this%dynconst
    CALL this%Info("            viscosity model:   " // this%viscosity%GetName())
    SELECT CASE(this%viscosity%GetType())
    CASE(molecular)
       WRITE (bulkconst_str,'(ES9.2)') this%bulkconst
       CALL this%Info("          dynamic viscosity:  " // trim(dynconst_str))
       CALL this%Info("             bulk viscosity:  " // trim(bulkconst_str))
    CASE(alpha,alpha_alt)
       CALL this%Info("                      alpha:  " // trim(dynconst_str))
    CASE(beta)
       CALL this%Info("                       beta:  " // trim(dynconst_str))
    CASE(powerlaw)
       WRITE (bulkconst_str,'(ES9.2)') this%power
       CALL this%Info("            coupling constant: " // trim(dynconst_str))
       CALL this%Info("                     exponent: " // trim(bulkconst_str))
    END SELECT
  END SUBROUTINE infosources


  SUBROUTINE updateviscosity(this,Mesh,Physics,Fluxes,time,pvar,cvar)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm, statevector_eulerisotherm
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Sources_viscosity),INTENT(INOUT) :: this
    CLASS(Mesh_base),INTENT(IN)         :: Mesh
    CLASS(Physics_base),INTENT(INOUT)   :: Physics
    CLASS(Fluxes_base),INTENT(IN)       :: Fluxes
    REAL,INTENT(IN)                     :: time
    CLASS(marray_compound),INTENT(INOUT)   :: pvar,cvar
    !------------------------------------------------------------------------!
    INTEGER           :: l
    !------------------------------------------------------------------------!
    ! only update viscosity if time has changed
    IF ((time.NE.this%time) .OR. (time .EQ. 0.)) THEN
       ! Set the dynamic viscosity
       SELECT CASE(this%viscosity%GetType())
       CASE(molecular)
          ! do nothing, since it has already been initialized
          ! in InitSources_viscosity
       CASE(alpha)
          ! compute Shakura-Sunyaev type alpha viscosity
          ! account for non-inertial frames by adding the frame angular velocity
          SELECT TYPE(p => pvar)
          CLASS IS(statevector_eulerisotherm)
            SELECT TYPE(phys => physics)
            CLASS IS(physics_eulerisotherm)
              this%dynvis%data1d(:) = etafkt_alpha(this%dynconst,p%density%data1d(:),&
                phys%bccsound%data1d(:),omega(this%ephir,p%velocity)+mesh%Omega)
            END SELECT
          END SELECT
       CASE(beta)
          ! compute Duschl type beta viscosity
          ! account for non-inertial frames by adding the frame angular velocity
          SELECT TYPE(p => pvar)
          CLASS IS(statevector_eulerisotherm)
            this%dynvis%data1d(:) = etafkt_beta(this%dynconst,mesh%radius%data2d(:,2), &
              p%density%data1d(:),omega(this%ephir,p%velocity)+mesh%Omega)
          END SELECT
       CASE(powerlaw)
          ! kinematic viscosity scales with power of specific angular momentum
          ! account for non-inertial frames by adding the frame angular velocity
          SELECT TYPE(p => pvar)
          CLASS IS(statevector_eulerisotherm)
            this%dynvis%data1d(:) = etafkt_powerlaw(this%dynconst,this%power, &
              mesh%radius%data2d(:,2),p%density%data1d(:),omega(this%ephir,p%velocity)+mesh%Omega)
          END SELECT
! TODO: height of the disk is calculated in another source module, which is not translated yet. Therefore, ALPHA_ALT is not
! working
!        CASE(ALPHA_ALT)
!          ! eta = nu*rho = alpha*cs*h * rho
!          ! compute dynamic viscosity
!          this%dynvis(:,:,:) = this%dynconst * Physics%bccsound%data3d(:,:,:) &
!                * Physics%sources%height(:,:,:) * pvar(:,:,:,Physics%DENSITY)
       END SELECT

       ! Set the bulk viscosity
       SELECT CASE(this%viscosity%GetType())
       CASE(molecular)
         ! do nothing, since it has already been initialized
         ! in InitSources_viscosity
       CASE(alpha,beta,powerlaw,alpha_alt)
         ! All available parametrized viscosities should be trace free.
         ! This can be achieved by setting the following:
         this%bulkvis%data1d(:) = -2./3.*this%dynvis%data1d(:)
       END SELECT

       ! update time
       this%time = time
    END IF
  CONTAINS
    ! some elemental functions for computation of the viscosity

    ! compute angular velocity
    PURE FUNCTION omega(ephir,velocity)
      IMPLICIT NONE
      !----------------------------------------------------------------------!
      CLASS(marray_base), INTENT(IN) :: ephir,velocity
      REAL, DIMENSION(SIZE(ephir%data2d,DIM=1)) :: omega
      !----------------------------------------------------------------------!
      INTEGER :: l
      !----------------------------------------------------------------------!
      ! project velocity on ephi/r -> local angular velocity
      ! use dynvis as temporary storage
      omega(:) = ephir%data2d(:,1) * velocity%data2d(:,1)
      DO l=2,ephir%DIMS(1)
        omega(:) = omega(:) + ephir%data2d(:,l) * velocity%data2d(:,l)
      END DO
    END FUNCTION omega

    ! Alpha viscosity
    ELEMENTAL FUNCTION etafkt_alpha(alpha,density,cs,omega) RESULT(eta)
      IMPLICIT NONE
      !----------------------------------------------------------------------!
      REAL, INTENT(IN) :: alpha,density,cs,omega
      REAL :: eta
      !----------------------------------------------------------------------!
      eta = alpha * density * cs**2 / abs(omega)
    END FUNCTION etafkt_alpha

    ! Beta viscosity
    ELEMENTAL FUNCTION etafkt_beta(beta,r,density,omega) RESULT(eta)
      IMPLICIT NONE
      !----------------------------------------------------------------------!
      REAL, INTENT(IN) :: beta,r,density,omega
      REAL :: eta
      !----------------------------------------------------------------------!
      ! nu  = beta * r * v_phi = beta * r^2 * omega
      ! eta = rho * nu
      eta = beta * density * r**2 * abs(omega)
    END FUNCTION etafkt_beta

    ! power law viscosity
    ELEMENTAL FUNCTION etafkt_powerlaw(beta,q,r,density,omega) RESULT(eta)
      IMPLICIT NONE
      !----------------------------------------------------------------------!
      REAL, INTENT(IN) :: beta,q,r,density,omega
      REAL :: eta
      !----------------------------------------------------------------------!
      eta = beta * density * (r**2 * abs(omega))**q
    END FUNCTION etafkt_powerlaw

  END SUBROUTINE updateviscosity


  SUBROUTINE externalsources_single(this,Mesh,Physics,Fluxes,Sources,time,dt,pvar,cvar,sterm)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm
    USE physics_euler_mod, ONLY : physics_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(sources_viscosity),INTENT(INOUT) :: this
    CLASS(mesh_base),INTENT(IN)            :: Mesh
    CLASS(physics_base),INTENT(INOUT)      :: Physics
    CLASS(fluxes_base),INTENT(IN)          :: Fluxes
    CLASS(sources_base), INTENT(INOUT)  :: Sources
    REAL,INTENT(IN)                        :: time, dt
    CLASS(marray_compound),INTENT(INOUT)   :: pvar,cvar,sterm
    !------------------------------------------------------------------------!
    CALL this%UpdateViscosity(mesh,physics,fluxes,time,pvar,cvar)
    SELECT TYPE(phys => physics)
    CLASS IS(physics_eulerisotherm)
      CALL phys%CalcStresses(mesh,pvar,this%dynvis,this%bulkvis, &
                this%btxx,this%btxy,this%btxz,this%btyy,this%btyz,this%btzz)
      CALL phys%ViscositySources(mesh,pvar,this%btxx,this%btxy,this%btxz, &
                this%btyy,this%btyz,this%btzz,sterm)
    END SELECT
  END SUBROUTINE externalsources_single


  SUBROUTINE calctimestep_single(this,Mesh,Physics,Fluxes,pvar,cvar,time,dt)
    USE physics_eulerisotherm_mod, ONLY : physics_eulerisotherm, statevector_eulerisotherm
    USE physics_euler_mod, ONLY : physics_euler, statevector_euler
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Sources_viscosity),INTENT(INOUT) :: this
    CLASS(Mesh_base),INTENT(IN)         :: Mesh
    CLASS(Physics_base),INTENT(INOUT)   :: Physics
    CLASS(Fluxes_base),INTENT(IN)       :: Fluxes
    CLASS(marray_compound), INTENT(INOUT) :: pvar,cvar
    REAL, INTENT(IN)                      :: time
    REAL, INTENT(OUT)                     :: dt
    !------------------------------------------------------------------------!
    REAL              :: invdt_x, invdt_y,invdt_z
    !------------------------------------------------------------------------!
    ! update dynamic and bulk viscosity
    CALL this%UpdateViscosity(mesh,physics,fluxes,time,pvar,cvar)

    ! compute kinematic viscosity
    SELECT TYPE(p => pvar)
    CLASS IS(statevector_eulerisotherm)
      this%kinvis%data1d(:) = this%dynvis%data1d(:) / p%density%data1d(:)
    END SELECT

    ! x-direction
    IF (mesh%INUM.GT.1) THEN
       invdt_x = maxval(this%kinvis%data1d(:) / mesh%dlx%data1d(:)**2, &
                        mask=mesh%without_ghost_zones%mask1d(:))
    ELSE
       ! set to zero, i.e. no limit in x-direction due to diffusion
       invdt_x = 0.0
    END IF
    ! y-direction
    IF (mesh%JNUM.GT.1) THEN
       invdt_y = maxval(this%kinvis%data1d(:)  / mesh%dly%data1d(:)**2, &
                        mask=mesh%without_ghost_zones%mask1d(:))
    ELSE
       ! set to zero, i.e. no limit in y-direction due to diffusion
       invdt_y = 0.0
    END IF
    ! z-direction
    IF (mesh%KNUM.GT.1) THEN
       invdt_z = maxval(this%kinvis%data1d(:)  / mesh%dlz%data1d(:)**2, &
                        mask=mesh%without_ghost_zones%mask1d(:))
    ELSE
       ! set to zero, i.e. no limit in z-direction due to diffusion
       invdt_z = 0.0
    END IF
    ! largest time step due to diffusion
    invdt_x = max(invdt_x, invdt_y,invdt_z)
    IF (invdt_x.GT.tiny(invdt_x)) THEN
       dt = this%cvis / invdt_x
    ELSE
       dt = huge(dt)
    END IF
  END SUBROUTINE calctimestep_single


  SUBROUTINE finalize(this)
    IMPLICIT NONE
    !------------------------------------------------------------------------!
    CLASS(Sources_viscosity),INTENT(INOUT) :: this
    !------------------------------------------------------------------------!
    IF (ALLOCATED(this%ephir)) THEN
      CALL this%ephir%Destroy()
      DEALLOCATE(this%ephir)
    END IF
    CALL this%dynvis%Destroy()
    CALL this%kinvis%Destroy()
    CALL this%bulkvis%Destroy()
    DEALLOCATE(this%dynvis,this%kinvis,this%bulkvis, &
               this%btxx,this%btyy,this%btzz,this%btxy,this%btxz,this%btyz)
  END SUBROUTINE finalize


END MODULE sources_viscosity_mod