! ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
MODULE WENOguts
CONTAINS
SUBROUTINE SmoothWeightPolyWENO(fs, isk, omega, qk, alpha, tmpSUM, fx)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
USE Interp_coeffs
! ----
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: fs(nvars,-enoOrder+1:enoOrder)
DOUBLE PRECISION, INTENT(OUT) :: fx(nvars)
DOUBLE PRECISION, INTENT(INOUT) :: isk(nvars,0:enoOrder)
DOUBLE PRECISION, INTENT(INOUT) :: qk(nvars,0:enoOrder)
DOUBLE PRECISION, INTENT(INOUT) :: omega(nvars,0:enoOrder)
DOUBLE PRECISION, INTENT(INOUT) :: tmpSUM(0:enoOrder,enoOrder)
DOUBLE PRECISION, INTENT(INOUT) :: alpha(0:enoOrder)
! ---- the regularization parameter
DOUBLE PRECISION, PARAMETER:: eps=1.D-06
! ---- WENO weighting parameter
INTEGER, PARAMETER:: pwr=2
INTEGER :: j,m,k
isk = 0.D0
DO j = 1,nvars,1
! ---- Compute inner SUM for m=1,enoOrder
tmpSUM = 0.D0
DO m = 1,enoOrder-1,1
DO k = 0,enoOrder,1
! ---- the coefficients dw
! ---- are held in the module Interp_coeffs
tmpSUM(k,m)=(SUM(dw(k,m,0:enoOrder-1)* &
& fs(j,-enoOrder+1+k:k)))**2
END DO
END DO
DO k = 0,enoOrder,1
isk(j,k) = SUM(tmpSUM(k,1:enoOrder-1))
END DO
! ---- make sure we don't miss our point with last stentcil.
isk(j,enoOrder) = maxval(isk(j,0:enoOrder))
END DO
DO j=1,nvars,1
! ---- use optimal weights (cw) and smoothness (isk)
! ---- to contruct the WENO reconstruction weights
alpha(0:enoOrder) = cw(0:enoOrder)/ &
& (eps+isk(j,0:enoOrder))**pwr
omega(j,0:enoOrder)= &
& alpha(0:enoOrder)/SUM(alpha(0:enoOrder))
END DO
DO k=0,enoOrder
DO j=1,nvars,1
qk(j,k)=SUM(aw(k,0:enoOrder-1)*fs(j,-enoOrder+1+k:k))
END DO
END DO
! ---- Convex superpostition of canidates
! ---- to make interpolates
DO j = 1, nvars, 1
fx(j) = sum( qk(j,0:enoOrder) * omega(j,0:enoOrder))
END DO
RETURN
END SUBROUTINE SmoothWeightPolyWENO
SUBROUTINE SignedFluxes(fps,fms,fs,us,alpha)
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ----
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: alpha(nvars)
DOUBLE PRECISION, INTENT(IN) :: fs(nvars,-enoOrder+1:enoOrder)
DOUBLE PRECISION, INTENT(IN) :: us(nvars,-enoOrder+1:enoOrder)
DOUBLE PRECISION, INTENT(OUT) :: fps(nvars,-enoOrder+1:enoOrder)
DOUBLE PRECISION, INTENT(OUT) :: fms(nvars,-enoOrder+1:enoOrder)
INTEGER :: l
DO l=-enoOrder+1, enoOrder
fps(:,l) = 0.5*(fs(:,l) + alpha(:)*us(:,l))
END DO
! ---- NOTICE: this loads backwards!
DO l=-enoOrder+1, enoOrder
fms(:,l) = 0.5*(fs(:,-l+1) - alpha(:)*us(:,-l+1))
END DO
END SUBROUTINE SignedFluxes
SUBROUTINE CalculateAlpha(alpha,lamda,lami,eta)
! ---- Alpha is used in the flux spliting.
! ---- using lami, eta are non-zero when
! ---- the Carbuncle fix is on.
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ----
IMPLICIT NONE
INTEGER :: i,j
DOUBLE PRECISION, INTENT(OUT) :: alpha(nvars)
DOUBLE PRECISION, INTENT(IN) :: lami(2,3)
DOUBLE PRECISION, INTENT(IN) :: eta
DOUBLE PRECISION, INTENT(IN) :: lamda(nvars)
DOUBLE PRECISION :: lamw(2,nvars)
lamw(:,1)=lami(:,3) ! u-c
do i=2,nvars-1
lamw(:,i)=lami(:,1) ! u
enddo
lamw(:,nvars)=lami(:,2) ! u+c
lamw=abs(lamw)
DO j=1,nvars,1
alpha(j)=max(dabs(lamda(j)), &
& maxval(lamw(:,j)),eta)
END DO
END SUBROUTINE CalculateAlpha
END MODULE WENOguts
! ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
MODULE OneDweno
CONTAINS
SUBROUTINE FluxInterpWENO(fxi,fx,ux,vx,lami,eta,dcflag,ilo,ihi,nn)
! ---- Subroutine uses the (physically) local values of a flux vector
! ---- fx defined over cells (-enoOrder+1:enoOrder) to calculate
! ---- the (WENO) interpolated flux, fxiW between the 0 and 1 cells.
! ---- Shared Variables
USE mesh
USE array_bounds
USE method_parms
! ---- Shared Procedures
USE Generic_EigenSystem
USE WENOguts
IMPLICIT NONE
include 'cles.i'
INTEGER, INTENT(IN) :: ilo, ihi, nn
DOUBLE PRECISION, INTENT(OUT) :: fxi(ncomps,ilo:ihi)
DOUBLE PRECISION, INTENT(IN) :: fx(nvars,ilo:ihi)
DOUBLE PRECISION, INTENT(IN) :: ux(ncomps,ilo:ihi)
DOUBLE PRECISION, INTENT(IN) :: vx(nvars,ilo:ihi)
INTEGER, INTENT(IN) :: dcflag(1:nn+1)
! ---- Carbuncle fix
DOUBLE PRECISION, INTENT(IN) :: eta(*)
DOUBLE PRECISION, INTENT(IN) :: lami(2,3,*)
! ---- local values projected into characterisitc space
DOUBLE PRECISION :: fs(nvars,-enoOrder+1:enoOrder)
DOUBLE PRECISION :: us(nvars,-enoOrder+1:enoOrder)
! ---- Plus (p) and Minus (m) going projected flux
DOUBLE PRECISION:: fps(nvars,-enoOrder+1:enoOrder)
DOUBLE PRECISION:: fms(nvars,-enoOrder+1:enoOrder)
! ---- Plus (p) and Minus (m) going interpolated flux
DOUBLE PRECISION :: fxp(nvars)
DOUBLE PRECISION :: fxm(nvars)
! ---- Eigen System
DOUBLE PRECISION :: lc(nvars,nvars)
DOUBLE PRECISION :: lamda(nvars)
DOUBLE PRECISION :: rc(nvars,nvars)
! ---- Smoothness measure
DOUBLE PRECISION :: isk(nvars,0:enoOrder)
! ---- WENO weights
DOUBLE PRECISION :: omega(nvars,0:enoOrder)
! ---- canidate interpolant values
DOUBLE PRECISION :: qk(nvars,0:enoOrder)
DOUBLE PRECISION :: tmpSUM(0:enoOrder,enoOrder)
DOUBLE PRECISION :: alphaW(0:enoOrder)
DOUBLE PRECISION:: alpha(nvars)
INTEGER :: i, l, k
INTEGER :: i0
! ---- ---- ---- ---- ---- ---- ---- ---- ----
! ---- Flux vector projection and split
! ---- ---- ---- ---- ---- ---- ---- ---- ----
do i=1, nn+1
if ( dcflag(i) .eq. CLES_SWITCH_WENO .or. use_dcflag .eq. CLES_FALSE ) then
i0 = i-1
! ---- Factor density-averaged flux matrix
CALL EigenSystem(ux(:,i0), ux(:,i), vx(:,i0), vx(:,i),lc, rc, lamda)
! ---- Project the local flux and vel onto the characteristics
DO k=-enoOrder+1,enoOrder,1
DO l=1,nvars,1
fs(l,k) = SUM(lc(l,:)*fx(:,i0+k))
us(l,k) = SUM(lc(l,1:nvars)*ux(1:nvars,i0+k))
END DO
END DO
if ( use_carbfix .eq. 1 ) then
! ---- Employ Carbuncle fix in finding alpha
CALL CalculateAlpha(alpha,lamda,lami(1,1,i),eta(i))
! ---- split flux into plus and minus directions
CALL SignedFluxes(fps,fms,fs,us,alpha)
else
DO l=-enoOrder+1, enoOrder
fps(1:nvars,l) = 0.5*(fs(1:nvars,l))
END DO
! ---- NOTICE: this loads backwards!
DO l=-enoOrder+1, enoOrder
fms(1:nvars,l) = 0.5*(fs(1:nvars,-l+1))
END DO
endif
if(stencil.eq.3) then
! just doing a simple flux splitting
! ---- Employ Carbuncle fix in finding alpha
CALL CalculateAlpha(alpha,lamda,lami(1,1,i),eta(i))
! ---- split flux into plus and minus directions
CALL SignedFluxes(fps,fms,fs,us,alpha)
fxp(1:nvars) = fps(1:nvars,0)
fxm(1:nvars) = fms(1:nvars,0)
else
! ---- ---- ---- ---- ---- ---- ---- ---- ----
! ---- Flux 'minus' interpolation
! ---- ---- ---- ---- ---- ---- ---- ---- ----
CALL SmoothWeightPolyWENO(fms,isk, omega, qk, alphaW, tmpSUM, fxm)
! ---- ---- ---- ---- ---- ---- ---- ---- ----
! ---- Flux 'plus' interpolation
! ---- ---- ---- ---- ---- ---- ---- ---- ----
CALL SmoothWeightPolyWENO(fps,isk, omega, qk, alphaW, tmpSUM, fxp)
! ---- Transform back to physical co-ordinates
end if
DO l=1,nvars,1
fxi(l,i) = SUM(rc(l,:)*(fxp(:)+fxm(:)))
END DO
endif
enddo
END SUBROUTINE FluxInterpWENO
END MODULE OneDweno