c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c Curvature based detection
c
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars,
$ ixlo,ixhi,iylo,iyhi, nx, ny, dx, dy,
$ direction,extra)
implicit none
include 'cles.i'
include 'cuser.i'
integer ncomps, nvars
integer ixlo, ixhi, iylo, iyhi
integer direction, extra(*)
INTEGER nx,ny
integer dcflag(1:nx+1,1:ny+1,2)
double precision ux(ncomps,ixlo:ixhi,iylo:iyhi)
double precision vx(nvars,ixlo:ixhi,iylo:iyhi)
DOUBLE PRECISION dx,dy
integer i, j, m, slow, shigh
double precision crvP, crvRho, crvS, x, y
double precision ThressCrvP,ThressCrvRho,ThressCrvS
integer dspan
dspan = 4
! ---- test criteria at each point
if ( direction .eq. 1 ) then
ThressCrvP = CurvP *dx*dx * 0.25d0
ThressCrvRho = CurvRho *dx*dx * 0.25d0
ThressCrvS = CurvS *dx*dx * 0.25d0
do j = 1, ny, 1
do i = 1, nx+1, 1
crvP = (vx(5,i-1,j)-2.*vx(5,i,j)+vx(5,i+1,j))
crvP = crvP/(vx(5,i-1,j)+2.*vx(5,i,j)+vx(5,i+1,j))
crvRho = (vx(1,i-1,j)-2.*vx(1,i,j)+vx(1,i+1,j))
crvRho = crvRho/(vx(1,i-1,j)+2.*vx(1,i,j)+vx(1,i+1,j))
crvS = (vx(6,i+1,j)-vx(6,i-1,j))/2.0d0
IF (((abs(crvP).gt. ThressCrvP).and.
$ (abs(crvrho).gt.ThressCrvRho)
$ .and.(crvP*crvRho.gt.0.0))) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(i-dspan,1)
shigh = min(i+dspan,nx+1)
do m=slow,shigh
dcflag(m,j,direction) = CLES_SWITCH_WENO
enddo
END IF
IF (abs(crvS).gt.CurvS) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(i-dspan,1)
shigh = min(i+dspan,nx+1)
do m=slow,shigh
dcflag(m,j,direction) = CLES_SWITCH_WENO
enddo
END IF
! Flag the boundaries
call cles_xLocation(i,x)
if ( x.lt.(GeomMin(1)+2.0d0*dx) ) then
if ( trunc .ne. 1 ) then
do m=1, i+1
dcflag(m,j,direction) = CLES_SWITCH_WENO
enddo
endif
else if (x.gt.(GeomMax(1)-2.0d0*dx)) then
do m=i,nx+1
dcflag(m,j,direction) = CLES_SWITCH_WENO
enddo
endif
enddo
enddo
else
ThressCrvP = CurvP *dy*dy * 0.25d0
ThressCrvRho = CurvRho *dy*dy * 0.25d0
ThressCrvS = CurvS *dy*dy * 0.25d0
do j = 1, ny+1, 1
DO i = 1, nx, 1
crvP = (vx(5,i,j-1)-2.*vx(5,i,j)+vx(5,i,j+1))
crvP = crvP/(vx(5,i,j-1)+2.*vx(5,i,j)+vx(5,i,j+1))
crvRho = (vx(1,i,j-1)-2.*vx(1,i,j)+vx(1,i,j+1))
crvRho = crvRho/(vx(1,i,j-1)+2.*vx(1,i,j)+vx(1,i,j+1))
crvS = (vx(6,i,j+1)-vx(6,i,j-1))/2.0d0
IF (((abs(crvP).gt. ThressCrvP).and.
$ (abs(crvrho).gt.ThressCrvRho)
$ .and.(crvP*crvRho.gt.0.0))) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(j-dspan,1)
shigh = min(j+dspan,ny+1)
do m=slow,shigh
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
END IF
IF (abs(crvS).gt.CurvS) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(j-dspan,1)
shigh = min(j+dspan,ny+1)
do m=slow,shigh
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
END IF
! Flag the boundaries
call cles_yLocation(j,y)
if ( y.lt.(GeomMin(2)+2.0d0*dy) ) then
do m=1, j+1
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
else if (y.gt.(GeomMax(2)-2.0d0*dy)) then
do m=j,ny+1
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
endif
enddo
enddo
endif
return
end