MODULE Generic_Transport
SAVE
! type of transport model
INTEGER :: trsp_visc_mode = 0
INTEGER :: trsp_kappa_mode = 0
INTEGER :: trsp_diff_mode = 0
! ---- the most viscous value in the flow
DOUBLE PRECISION :: maxMu = 0.0d0
DOUBLE PRECISION :: maxKappa = 0.0d0
DOUBLE PRECISION :: maxDiff = 0.0d0
! this is esentially (gamma-1)/gamma for air used in computing
! maxKappa from the pressure and temperature. This is needed to
! evaluate Kappa/rho Cp (approximately).
DOUBLE PRECISION, PRIVATE :: kfactor
PARAMETER (kfactor = 0.4d0/1.4d0)
! ---- used for the viscous terms
DOUBLE PRECISION, ALLOCATABLE :: mu(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: kappa(:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: diff(:,:,:,:)
! --- These routines maybe altered if we want different
! --- base viscosities for the different gasses
INTERFACE GetTransport
MODULE PROCEDURE GetTransportOneD
MODULE PROCEDURE GetTransportTwoD
MODULE PROCEDURE GetTransportThreeD
END INTERFACE
interface
subroutine cles_transport(ux, ncomps, vx, nvars, visc, kappa, &
diff, nscal, n, ierr)
integer nvars, ncomps, n, nscal, ierr
double precision ux(ncomps,*), vx(nvars,*)
double precision visc(*), kappa(*), diff(*)
end subroutine cles_transport
end interface
private cles_transport
CONTAINS
SUBROUTINE SetupTransport()
use method_parms
IMPLICIT NONE
include 'cles.i'
integer ierr
interface
subroutine cles_setup_transport(visc_mode, kappa_mode, diff_mode, ierr)
integer visc_mode, kappa_mode, diff_mode, ierr
end subroutine cles_setup_transport
end interface
call cles_setup_transport(trsp_visc_mode, trsp_kappa_mode, &
trsp_diff_mode, ierr)
RETURN
END SUBROUTINE SetupTransport
SUBROUTINE GetTransportOneD(ux,vx)
! ---- Shared Variables
USE mesh
use array_bounds
USE method_parms
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi)
INTEGER n, ierr
n = ixhi-ixlo+1
call cles_transport(ux, ncomps, vx, nvars, mu, kappa, diff, nscal, n, ierr)
! --- this is used when computing the timestep size
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
maxMu = maxval(mu(1:nx,1,1)/vx(1,1:nx))
else
maxMu = mu(1,1,1)/minval(vx(1,1:nx))
endif
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
maxKappa = maxval(kappa(1:nx,1,1)*ux(nvars+1,1:nx)/vx(5,1:nx))
else
maxKappa = kappa(1,1,1)*maxval(ux(nvars+1,1:nx)/vx(5,1:nx))
endif
maxKappa = maxKappa*kfactor
END SUBROUTINE GetTransportOneD
SUBROUTINE GetTransportTwoD(ux,vx)
! ---- Shared Variables
USE mesh
use array_bounds
USE method_parms
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi)
INTEGER n, ierr
n = (ixhi-ixlo+1)*(iyhi-iylo+1)
call cles_transport(ux, ncomps, vx, nvars, mu, kappa, diff, nscal, n, ierr)
! --- this is used when computing the timestep size
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
maxMu = maxval(mu(1:nx,1:ny,1)/vx(1,1:nx,1:ny))
else
maxMu = mu(1,1,1)/minval(vx(1,1:nx,1:ny))
endif
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
maxKappa = maxval(kappa(1:nx,1:ny,1)*ux(nvars+1,1:nx,1:ny) &
/vx(5,1:nx,1:ny))
else
maxKappa = kappa(1,1,1)*maxval(ux(nvars+1,1:nx,1:ny) &
/vx(5,1:nx,1:ny))
endif
maxKappa = maxKappa*kfactor
END SUBROUTINE GetTransportTwoD
SUBROUTINE GetTransportThreeD(ux,vx)
! ---- Shared Variables
USE mesh
use array_bounds
USE method_parms
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER n, ierr
call cleslog_log_enter('GetTransportThreeD')
n = (ixhi-ixlo+1)*(iyhi-iylo+1)*(izhi-izlo+1)
call cles_transport(ux, ncomps, vx, nvars, mu, kappa, diff, nscal, n, ierr)
! --- this is used when computing the timestep size
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
maxMu = maxval(mu(1:nx,1:ny,1:nz)/vx(1,1:nx,1:ny,1:nz))
else
maxMu = mu(1,1,1)/minval(vx(1,1:nx,1:ny,1:nz))
endif
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
maxKappa = maxval(kappa(1:nx,1:ny,1:nz)*ux(nvars+1, &
1:nx,1:ny,1:nz)/vx(5,1:nx,1:ny,1:nz))
else
maxKappa = kappa(1,1,1)*maxval(ux(nvars+1, &
1:nx,1:ny,1:nz) /vx(5,1:nx,1:ny,1:nz))
endif
maxKappa = maxKappa*kfactor
call cleslog_log_exit('GetTransportThreeD')
END SUBROUTINE GetTransportThreeD
SUBROUTINE AllocateTransport()
! ---- Shared variables used in dimensioning
USE mesh
use array_bounds
USE method_parms
IMPLICIT NONE
include 'cles.i'
integer iret
call cleslog_log_enter('AllocateTransport')
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
if ( dim .eq. 1 ) then
ALLOCATE ( mu(ixlo:ixhi,1,1), stat = iret )
else if ( dim .eq. 2 ) then
ALLOCATE ( mu(ixlo:ixhi,iylo:iyhi,1), stat = iret )
else
ALLOCATE ( mu(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat = iret )
endif
else
ALLOCATE ( mu(1,1,1), stat = iret )
endif
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AllocateTransport/mu', iret)
endif
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
if ( dim .eq. 1 ) then
ALLOCATE ( kappa(ixlo:ixhi,1,1), stat = iret )
else if ( dim .eq. 2 ) then
ALLOCATE ( kappa(ixlo:ixhi,iylo:iyhi,1), stat = iret )
else
ALLOCATE ( kappa(ixlo:ixhi,iylo:iyhi,izlo:izhi), stat = iret )
endif
else
ALLOCATE ( kappa(1,1,1), stat = iret )
endif
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AllocateTransport/kappe', iret)
endif
if ( nscal .gt. 0 ) then
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
if ( dim .eq. 1 ) then
ALLOCATE ( diff(ixlo:ixhi,nscal,1,1), stat = iret )
else if ( dim .eq. 2 ) then
ALLOCATE ( diff(ixlo:ixhi,iylo:iyhi,nscal,1), stat = iret )
else
ALLOCATE ( diff(ixlo:ixhi,iylo:iyhi,izlo:izhi,nscal), stat = iret )
endif
else
ALLOCATE ( diff(nscal,1,1,1), stat = iret)
endif
else
ALLOCATE ( diff(1,1,1,1), stat = iret)
endif
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'AllocateTransport/diff', iret)
endif
call cleslog_log_exit('AllocateTransport')
END SUBROUTINE AllocateTransport
SUBROUTINE DeallocateTransport()
! ---- Shared variables used in dimensioning
USE mesh
USE method_parms
IMPLICIT NONE
include 'cles.i'
integer iret
call cleslog_log_enter('DeallocateTransport')
DEALLOCATE ( mu, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'DeallocateTransport/mu', iret)
endif
DEALLOCATE ( kappa, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'DeallocateTransport/kappa', iret)
endif
DEALLOCATE ( diff, stat = iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'DeallocateTransport/diff', iret)
endif
call cleslog_log_exit('DeallocateTransport')
END SUBROUTINE DeallocateTransport
END Module Generic_Transport
SUBROUTINE Transport_Viscosity(var)
use mesh
use array_bounds
use method_parms
use Generic_Transport
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(OUT):: var(ixlo:ixhi,iylo:iyhi,izlo:izhi)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
var = mu
else
var = mu(1,1,1)
endif
END SUBROUTINE Transport_Viscosity
SUBROUTINE Transport_Kappa(var)
use mesh
use array_bounds
use method_parms
use Generic_Transport
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(OUT):: var(ixlo:ixhi,iylo:iyhi,izlo:izhi)
if ( trsp_kappa_mode .eq. CLES_TRANSP_KAPPA_VAR ) then
var = kappa
else
var = kappa(1,1,1)
endif
END SUBROUTINE Transport_Kappa
SUBROUTINE Transport_Diffusion(var, is)
use mesh
use array_bounds
use method_parms
use Generic_Transport
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: is
DOUBLE PRECISION, INTENT(OUT):: var(ixlo:ixhi,iylo:iyhi,izlo:izhi)
if ( nscal .gt. 0 ) then
if ( trsp_diff_mode .eq. CLES_TRANSP_DIFF_VAR ) then
var = diff(:,:,:,is)
else
var = diff(is,1,1,1)
endif
endif
END SUBROUTINE Transport_Diffusion