c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c Non-linear mapping/entropy condition based detection
c
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars,
$ ixlo,ixhi,nx,dx,direction,extra)
implicit none
include 'airsf6.i'
include 'cles.i'
include 'cuser.i'
integer ncomps, nvars
integer ixlo, ixhi
integer nx, direction, extra(*)
integer dcflag(1:nx+1,1)
double precision ux(ncomps,ixlo:ixhi)
double precision vx(nvars,ixlo:ixhi)
double precision dx
integer i, m, slow, shigh, l
double precision dl, x
double precision gl, Rol, ul, asl, hl, rhosqrtl
double precision gr, Ror, ur, asr, hr, rhosqrtr
double precision gtilde, utilde, astilde, htilde, rhotilde
double precision mu, tmp, scalg, scalm
double precision d1r, d1l, d3r, d3l, theta
integer dcspan, nfield
parameter (dcspan = 3)
parameter (nfield = 5) ! pressure
dl = dx
c--------------------------------------------------------------------------
c Roe's linearized Riemann problem based detection
do i = 1, nx+1, 1
c First, determine the states for the local linearized Riemann problem
c at the interface "i" between the cell centers "i-1" and "i".
c Left state ql=q(i-1)
call cles_gamma(vx(1,i-1),nvars,gl,Rol)
ul = vx(2,i-1)
asl = sqrt(gl*vx(5,i-1)/vx(1,i-1))
hl = (ux(5,i-1)+vx(5,i-1))/vx(1,i-1)
rhosqrtl = sqrt(vx(1,i-1))
c Right state qr=q(i)
call cles_gamma(vx(1,i),nvars,gr,Ror)
ur = vx(2,i)
asr = sqrt(gr*vx(5,i)/vx(1,i))
hr = (ux(5,i)+vx(5,i))/ux(1,i)
rhosqrtr = sqrt(vx(1,i))
c Roe-averaged quantities
call cles_roe(ux(1,i-1),ux(1,i),ncomps,vx(1,i-1),vx(1,i),
& nvars,gtilde,tmp,mu,0)
utilde = (rhosqrtl*vx(2,i-1) + rhosqrtr*vx(2,i))
& /(rhosqrtl + rhosqrtr)
htilde = (rhosqrtl*hl + rhosqrtr*hr)
& /(rhosqrtl + rhosqrtr)
astilde = sqrt((gtilde-1.0d0)*(htilde-0.5d0*utilde**2))
rhotilde = sqrt(vx(1,i-1)*vx(1,i))
c Second, apply normalized Lax's entropy condition (with some
c departure 'LaxThres' allowed) to determine which wave is a
c shock for the Riemann problem at the cell interface "i".
c Theoretically, a shock wave can only have the Jacobian eigenvalue
c u+a or u-a.
c Lax's residual for wave with speed u-a
d1r = (utilde-astilde)-(ur-asr)
d1r = d1r/astilde
d1l = (ul-asl)-(utilde-astilde)
d1l = d1l/astilde
c Lax's residual for wave with speed u+a
d3r = (utilde+astilde)-(ur+asr)
d3r = d3r/astilde
d3l = (ul+asl)-(utilde+astilde)
d3l = d3l/astilde
c preselect strength
if ((d1r.gt.0.0d0).and.(d1l.gt.0.0d0).and.
& (d3r.gt.0.0d0).and.(d3l.gt.0.0d0)) then !preselect shock
if (((d1r.gt.ShockLaxThres).and.(d1l.gt.ShockLaxThres))
& .or.((d3r.gt.ShockLaxThres).and.(d3l.gt.ShockLaxThres))) then
c Third, flag shock using a mapping of pressure profile,
c with some thresholds to avoid flagging very weak shocks
theta = abs((vx(nfield,i)-vx(nfield,i-1))
& /(vx(nfield,i)+vx(nfield,i-1)))
theta = 2.0d0*theta/(1.0d0+theta)**2.0d0
if (theta.gt.MappThres) then
slow = max(i-dcspan,1)
shigh = min(i+1+dcspan,nx+1)
do m = slow, shigh
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
endif
endif
endif
scalm = (vx(6,i) + vx(6,i-1))*0.5d0
scalg = abs(vx(6,i) - vx(6,i-1))
if ( scalg .gt. ScalThres .or. scalm .lt. -1.01d0
$ .or. scalm .gt. 1.01d0 ) then
slow = max(i-2*dcspan,1)
shigh = min(i+1+2*dcspan,nx+1)
do m = slow, shigh
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
endif
enddo
! Flag the wall
do i = 1, nx
call cles_xLocation(i, x)
if ( x .gt. xmax-3*dx ) then
dcflag(i+1,direction) = CLES_SWITCH_WENO
endif
enddo
return
end