c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c
c Curvature based detection
c
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine cles_dcflag(ux,vx,dcflag,ncomps,nvars,
$ ixlo,ixhi,nx,dx,direction,extra)
implicit none
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, nspan
double precision grdS1, grdS2, x
double precision apress, bpress, cpress
double precision arho, brho, crho
integer dcspan
parameter (dcspan=6)
c--------------------------------------------------------------------------
c Use gradient detection for the scalar
c test criteria at each point
do i = 1, nx, 1
grdS1 = (vx(6,i+1)-vx(6,i-1))/2.0d0
grdS2 = (vx(7,i+1)-vx(7,i-1))/2.0d0
if (abs(grdS1).gt.gradS .or. abs(grdS2).gt.gradS) then
c we will use WENO here
slow = max(i-dcspan,1)
shigh = min(i+1+dcspan,nx+1)
do m=slow,shigh
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
end if
enddo
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
c Correlation detection
c This is esentially (nghost - span)*2.
nspan = (1-ixlo-span)*2
if ( nspan .lt. 2 ) then
print *, 'Number of ghost cells is too low for'//
$ ' non-linear projection flagging'
stop
endif
do l=1-nspan, 1+nx*2+nspan, 1
if ( mod(l,2) .eq. 0 ) then
c Center of cell
i = l/2
c Correlation of the pressure field (nfield=5) with a centered tanh
call TanhCorrel(span,span,thvec,sumth,sumth2,
& vx,nvars,ixlo,ixhi,i,5,apress,bpress,cpress)
if ( abs(bpress) .lt. minjump*abs(apress) ) cycle
if ( abs(cpress).lt. minc) cycle
slow = max(i-dcspan,1)
shigh = min(i+dcspan+1,nx+1)
do m=slow,shigh
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
else
c Wall of cell
i = 1+l/2
c Correlation of the pressure field (nfield=5) with a shifted tanh
call TanhCorrel(span,span-1,thvecbis,sumthbis,sumth2bis,
& vx,nvars,ixlo,ixhi,i,5,apress,bpress,cpress)
if ( abs(bpress) .lt. minjump*abs(apress) ) cycle
if ( abs(cpress).lt. minc) cycle
slow = max(i-dcspan,1)
shigh = min(i+dcspan,nx+1)
do m=slow,shigh
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
endif
enddo
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
c Flag the boundaries
do i=1, nx
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,direction) = CLES_SWITCH_WENO
enddo
endif
else if (x.gt.(GeomMax(1)-2.0d0*dx)) then
do m=i,nx+1
dcflag(m,direction) = CLES_SWITCH_WENO
enddo
endif
enddo
c--------------------------------------------------------------------------
return
end
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine TanhCorrel(span,span1,thvec,sumth,sumth2,
& vx,nvars,ixlo,ixhi,i,nfield,a,b,c)
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
c Correlation based detection of the field vx(.,nfield) with a tanh
c profile. Local least-square fitting of a+b*tanh(x/dwidth) or abis
c +bbis*tanh((x+dx/2)/dwidth)
c Thickness of the profile in terms of grid sizes. A value of 1 means
c a discontinuity of 6 grids points because of the thickness of the
c tanh(): dwidth (e.g. = 0.75d0)
c Minimum value of a correlation coefficient considered to be a good
c match: minc (e.g. = 0.9d0)
c To discriminate small oscillations in the coorelated field that have
c no meaning, we stipulate that any jump in the field below 'minjump'
c of the mean value is not a jump (e.g. minjump = 0.02d0 for 2%)
implicit none
c include 'cuser.i'
double precision sumth,sumth2
integer span,span1,nvars,ixlo,ixhi,i,nfield
double precision thvec(-span:span1)
double precision vx(nvars,ixlo:ixhi)
double precision a,b,c
integer m
double precision sumy,sumyt,sum
double precision det,dp
a = 0.0d0
b = 0.0d0
det = 0.0d0
sumy = 0.0d0
sumyt = 0.0d0
do m=-span,span1, 1
sumy = sumy + vx(nfield,i+m)
sumyt = sumyt + vx(nfield,i+m)*thvec(m)
enddo
c Optimal values
det = dble(span1+span+1)*sumth2-sumth*sumth
a = (sumy*sumth2-sumth*sumyt)/det
b = (dble(span1+span+1)*sumyt-sumth*sumy)/det
c Compute correlation
dp = 0.0d0
c = 0.0d0
sum = 0.0d0
sumy = 0.0d0
do m=-span,span1, 1
dp = vx(nfield,i+m)-a
sum = sum + dp*thvec(m)
sumy = sumy + dp*dp
enddo
if ( abs(sumy) .lt. 1.0d-12 ) then
c = 0.0d0
else
c = sum/sqrt(sumy*sumth2)
endif
return
end