c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c Flux smoothness based detection
c
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars,
$ ixlo,ixhi,iylo,iyhi,izlo,izhi,nx,ny,nz,dx,dy,dz,
$ direction,extra)
implicit none
integer ncomps, nvars
integer ixlo, ixhi, iylo, iyhi, izlo, izhi
integer direction, extra(*)
INTEGER nx,ny,nz
integer dcflag(1:nx+1,1:ny+1,1:nz+1,3)
double precision ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
double precision vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION dx,dy,dz
integer i, j, k, m, slow, shigh, loop, size, enoOrder
double precision smoothness, dmdx, maxv(4), minv(4), variance(4)
integer span
double precision TriggerValue, smooth_eps
TriggerValue = 0.0d0
smooth_eps = 0.001
span = 3
enoOrder = extra(1)
size = enoOrder-1
! ---- test criteria at each point
if ( direction .eq. 1 ) then
do k = 1, nz, 1
do j = 1, ny, 1
DO i = 1, nx+1, 1
DO loop =1,4
do m=i-size, i+size-1
dmdx = abs(fx(loop,m+1,j,k)
$ -fx(loop,m-1,j,k))
if ( m .eq. i-size ) then
maxv(loop) = dmdx
minv(loop) = dmdx
else
maxv(loop) = MAX(maxv(loop), dmdx)
minv(loop) = MIN(minv(loop), dmdx)
endif
enddo
variance(loop) = maxv(loop)/(minv(loop)+smooth_eps)
if ( loop .eq. 1 ) then
smoothness = variance(loop)
else
smoothness = MAX(smoothness, variance(loop))
endif
END DO
IF ( smoothness .gt. TriggerValue ) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(i-span,1)
shigh = min(i+span,nx+1)
do m=slow,shigh
dcflag(m,j,k,direction) = 1
enddo
END IF
enddo
enddo
enddo
else if (direction .eq. 2 ) then
do k = 1, nz, 1
do j = 1, ny+1, 1
do i = 1, nx, 1
DO loop =1,4
do m=j-size, j+size-1
dmdx = abs(fx(loop,i,m+1,k)
$ -fx(loop,i,m-1,k))
if ( m .eq. j-size ) then
maxv(loop) = dmdx
minv(loop) = dmdx
else
maxv(loop) = MAX(maxv(loop), dmdx)
minv(loop) = MIN(minv(loop), dmdx)
endif
enddo
variance(loop) = maxv(loop)/(minv(loop)+smooth_eps)
if ( loop .eq. 1 ) then
smoothness = variance(loop)
else
smoothness = MAX(smoothness, variance(loop))
endif
END DO
IF ( smoothness .gt. TriggerValue ) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(j-span,1)
shigh = min(j+span,ny+1)
do m=slow,shigh
dcflag(i,m,k,direction) = 1
enddo
END IF
enddo
enddo
enddo
else
do k = 1, nz+1, 1
do j = 1, ny, 1
DO i = 1, nx, 1
DO loop =1,4
do m=k-size, k+size-1
dmdx = abs(fx(loop,i,j,m+1)
$ -fx(loop,i,j,m-1))
if ( m .eq. k-size ) then
maxv(loop) = dmdx
minv(loop) = dmdx
else
maxv(loop) = MAX(maxv(loop), dmdx)
minv(loop) = MIN(minv(loop), dmdx)
endif
enddo
variance(loop) = maxv(loop)/(minv(loop)+smooth_eps)
if ( loop .eq. 1 ) then
smoothness = variance(loop)
else
smoothness = MAX(smoothness, variance(loop))
endif
END DO
IF ( smoothness .gt. TriggerValue ) THEN
! we are at a shock
! ---- we will use WENO here
slow = max(k-span,1)
shigh = min(k+span,nz+1)
do m=slow,shigh
dcflag(i,j,m,direction) = 1
enddo
END IF
enddo
enddo
enddo
endif
return
end