SUBROUTINE SetInterpConstants
! ---- This subroutine sets the values of the interpolation
! ---- coeffs used in WENO and in Center Difference
! ---- It assumes that the 'stencil' has already been set.
! ---- 'stencil' is shared in Interp_coeffs
! ---- Also assumes that the arrays which hold the interpolation
! ---- coeffs have already been established.
! ---- Shared Variables
USE Interp_coeffs ! --- shares: cw,dw,aw,cd1,cd2...
USE method_parms ! --- shares: nghost, enoOrder,stencil..
! ----
IMPLICIT NONE
! ---- this is the same as /alpha in Hill and Pullin JCP
DOUBLE PRECISION :: bandwidth_par
IF ( stencil .lt. 3 .or. stencil .gt. 7 ) THEN
PRINT *, 'You have requested a ',stencil,' point stencil'
PRINT *, 'but this code only works for 3 or 7 point stencils'
STOP
END IF
! ---- initialize the WENO coeffs as zero
aw = 0.0d0
cw = 0.0d0
dw = 0.0d0
ac=1./6. ; cc=1./10. ; dc=1./2. ; b2=sqrt(13./12.)
! WENO Stencil part
if ( stencil .le. 5 ) then
! ---- coeffs for the (WENO) canidate interpolation stencils.
aw(0,0) = -0.5d0
aw(0,1) = 3.0d0/2.0d0
aw(1,0) = 0.5d0
aw(1,1) = 0.5d0
aw(2,0) = 3.0d0/2.0d0
aw(2,1) = -0.5d0
! ---- weights for combining canidates stencils: Standard WENO
cw(0) = 1.0d0/3.0d0
cw(1) = 2.0d0/3.0d0
cw(2) = 0.0d0
! The center difference stencil
dw(0,1,0) = 1.0d0
dw(0,1,1) = -1.0d0
dw(1,1,0) = 1.0d0
dw(1,1,1) = -1.0d0
dw(2,1,0) = 1.0d0
dw(2,1,1) = -1.0d0
else
! ---- coeffs for the (WENO) canidate interpolation stencils.
aw(0,0) = 2.*ac
aw(0,1) = -7.*ac
aw(0,2) = 11.*ac
aw(1,0) = -ac
aw(1,1) = 5.*ac
aw(1,2) = 2.*ac
aw(2,0) = 2.*ac
aw(2,1) = 5.*ac
aw(2,2) = -ac
! ---- additional canidate stencil: for optimized WENO
aw(3,0) = 11.*ac
aw(3,1) = -7.*ac
aw(3,2) = 2.*ac
! ---- weights for combining canidates stencils: Standard WENO
cw(0) = cc
cw(1) = 6.0*cc
cw(2) = 3.0*cc
cw(3) = 0.0
! ---- coeffs used in determining smoothness inside WENO
! ---- used in canidate weighting
dw(0,1,0) = dc
dw(0,1,1) = -4.*dc
dw(0,1,2) = 3.*dc
dw(1,1,0) = -dc
dw(1,1,1) = 0.
dw(1,1,2) = dc
dw(2,1,0) = 3.*dc
dw(2,1,1) = -4.*dc
dw(2,1,2) = dc
! ---- these 3 for optimized
dw(3,1,0) = -5.*dc
dw(3,1,1) = 8.*dc
dw(3,1,2) = -3.*dc
dw(0,2,0) = b2
dw(0,2,1) = -2.*b2
dw(0,2,2) = b2
dw(1,2,0) = b2
dw(1,2,1) = -2.*b2
dw(1,2,2) = b2
dw(2,2,0) = b2
dw(2,2,1) = -2.*b2
dw(2,2,2) = b2
! ---- these 3 for optimised
dw(3,2,0) = b2
dw(3,2,1) = -2.*b2
dw(3,2,2) = b2
endif
whichstencilsize: SELECT CASE(stencil)
CASE (3)
! ---- coeffs for center diff
tcd_center(3) = 0.0d0
tcd_center(2) = 0.0d0
tcd_center(1) = 0.5d0
tcd_bndry(1,1) = -1.0d0
tcd_bndry(2,1) = 1.0d0
tcd_bndry(3,1) = 0.0d0
tcd_bndry(4,1) = 0.0d0
tcd_bndry(1,2) = -0.5d0
tcd_bndry(2,2) = 0.0d0
tcd_bndry(3,2) = 0.5d0
tcd_bndry(4,2) = 0.0d0
! ---- weights for combining canidates stencils: optimized WENO
cw(0) = 0.0d0
cw(1) = 1.0d0
cw(2) = 0.0d0
tcd_flux(1) = tcd_center(1)
tcd_flux(2) = 0.0d0
tcd_flux(3) = 0.0d0
CASE (5)
! ---- five point stencil
! ---- This is the free parameter determined by Ghosal
! ---- Change this to get other stencils.
! ---- set to -1.d0/12.d0 for standard 4-th order stencil
if ( order .eq. 4 ) then
bandwidth_par = -1.0d0/12.0d0
else
! optimized
bandwidth_par = -0.197d0
endif
! ---- coeffs for center diff
tcd_center(3) = 0.0d0
tcd_center(2) = bandwidth_par
tcd_center(1) = 0.50d0 - 2.0d0*bandwidth_par
tcd_bndry(1,1) = -0.5d0/(0.50d0+tcd_center(2))
tcd_bndry(2,1) = (0.5d0-tcd_center(2))/(0.5d0+tcd_center(2))
tcd_bndry(3,1) = tcd_center(2)/(0.5d0+tcd_center(2))
tcd_bndry(4,1) = 0.0d0
tcd_bndry(1,2) = -(0.5d0-tcd_center(2))/(1.0d0-tcd_center(2))
tcd_bndry(2,2) = 0.0d0
tcd_bndry(3,2) = tcd_center(1)/(1.0d0-tcd_center(2))
tcd_bndry(4,2) = tcd_center(2)/(1.0d0-tcd_center(2))
! ---- weights for combining canidates stencils: optimized WENO
cw(0) = -2.0d0 * bandwidth_par
cw(1) = 1.0d0 + 2.0d0*bandwidth_par
cw(2) = cw(0)
tcd_flux(3) = 0.0d0
tcd_flux(2) = tcd_center(2)
tcd_flux(1) = tcd_flux(2)+tcd_center(1)
CASE (7)
! ---- seven point stencil
! ---- This is the free parameter determined by Ghosal
! ---- Change this to get other stencils.
! ---- set to 1.d0/60.d0 for standard 6-th order stencil
if ( order .eq. 6 ) then
bandwidth_par = 1.d0/60.d0
else
! optimized
bandwidth_par = 0.0605d0
endif
! ---- coeffs for center diff
tcd_center(3) = bandwidth_par
tcd_center(2) = - 1.d0/12.d0 - 4.0*bandwidth_par
tcd_center(1) = 2.d0/3.d0 + 5.*bandwidth_par
tcd_bndry(1,1) = -0.5d0/(0.50d0+tcd_center(2))
tcd_bndry(2,1) = (0.5d0-tcd_center(2))/(0.5d0+tcd_center(2))
tcd_bndry(3,1) = tcd_center(2)/(0.5d0+tcd_center(2))
tcd_bndry(4,1) = 0.0d0
tcd_bndry(1,2) = -(0.5d0-tcd_center(2))/(1.0d0-tcd_center(2))
tcd_bndry(2,2) = 0.0d0
tcd_bndry(3,2) = tcd_center(1)/(1.0d0-tcd_center(2))
tcd_bndry(4,2) = tcd_center(2)/(1.0d0-tcd_center(2))
! ---- weights for combining canidates stencils: optimized WENO
cw(0) = 3.0*tcd_center(3)
cw(1) = 0.5d0 - cw(0)
cw(2) = cw(1)
cw(3) = cw(0)
tcd_flux(3) = tcd_center(3)
tcd_flux(2) = tcd_flux(3)+tcd_center(2)
tcd_flux(1) = tcd_flux(2)+tcd_center(1)
END SELECT whichstencilsize
END SUBROUTINE SetInterpConstants
SUBROUTINE TestGhostSize
USE method_parms ! --- shares: nghost, enoOrder,stencil,Thresshold..
IMPLICIT NONE
include 'cles.i'
! ---- The number of ghost cells
! ---- For LES or viscous: nghost >= enoOrder + 1
! ---- otherwise nghost >= enoOrder
IF (useViscous.eq.CLES_TRUE) THEN
IF(nghost.lt.enoOrder + 1) THEN
PRINT *, 'Set GhostCells (in the GridHierarchy{}-section) to ',enoOrder+1,'!'
STOP
END IF
ELSE
IF(nghost.lt.enoOrder) THEN
PRINT *, 'Set GhostCells (in the GridHierarchy{}-section) to ',enoOrder+1,'!'
STOP
END IF
END IF
END SUBROUTINE TestGhostSize
SUBROUTINE SetRKConstants(rk, c_rk1, c_rk2, c_rk3,f_rk)
IMPLICIT NONE
INTEGER, INTENT(IN) :: rk
DOUBLE PRECISION, INTENT(OUT) :: c_rk1, c_rk2, c_rk3,f_rk
call SetRKConstants_SSP(rk, c_rk1, c_rk2, c_rk3,f_rk)
RETURN
END SUBROUTINE SetRKConstants
SUBROUTINE SetRKConstants_SSP(rk, c_rk1, c_rk2, c_rk3,f_rk)
IMPLICIT NONE
INTEGER, INTENT(IN) :: rk
DOUBLE PRECISION, INTENT(OUT) :: c_rk1, c_rk2, c_rk3,f_rk
! --- for 3rd order SSP Runge-Kutta
whichRK: SELECT CASE (rk)
CASE (1)
c_rk1 = 1.0d0
c_rk2 = 0.0d0
c_rk3 = 1.0d0
f_rk = 1.d0/6.0d0
CASE (2)
c_rk1 = 0.25d0
c_rk2 = 0.25d0
c_rk3 = 0.75d0
f_rk = 1.d0/6.0d0
CASE (3)
c_rk1 = 2.0d0/3.0d0
c_rk2 = 2.0d0/3.0d0
c_rk3 = 1.0d0/3.0d0
f_rk = 4.d0/6.0d0
END SELECT WHICHRK
RETURN
END SUBROUTINE SetRKConstants_SSP