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