! ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 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