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 nx, ny, direction, extra(*)
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, l, nspanx, nspany
double precision grdS1, grdS2, x, y
double precision apress, bpress, cpress
double precision arho, brho, crho
integer dcspan
parameter (dcspan=6)
c--------------------------------------------------------------------------
c Use gradient detection for the scalar
if ( direction .eq. 1 ) then
c test criteria at each point
do j=1, ny, 1
do i = 1, nx, 1
grdS1 = (vx(6,i+1,j)-vx(6,i-1,j))/2.0d0
grdS2 = (vx(7,i+1,j)-vx(7,i-1,j))/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,j,direction) = CLES_SWITCH_WENO
enddo
end if
enddo
enddo
else
do j=1, ny, 1
do i = 1, nx, 1
grdS1 = (vx(6,i,j+1)-vx(6,i,j-1))/2.0d0
grdS2 = (vx(7,i,j+1)-vx(7,i,j-1))/2.0d0
if (abs(grdS1).gt.gradS .or.
$ abs(grdS2).gt.gradS) then
c we will use WENO here
slow = max(j-dcspan,1)
shigh = min(j+1+dcspan,ny+1)
do m=slow,shigh
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
end if
enddo
enddo
endif
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
c Correlation detection
if ( direction .eq. 1 ) then
c This is esentially (nghost - span)*2.
nspanx = (1-ixlo-span)*2
if ( nspanx .lt. 2 ) then
print *, 'Number of ghost cells is too low for'//
$ ' non-linear projection flagging'
stop
endif
do j=1, ny, 1
do l=1-nspanx, 1+nx*2+nspanx, 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,iylo,iyhi,i,j,5,direction,
$ apress,bpress,cpress)
if ( abs(bpress) .lt. minjump*abs(apress) ) cycle
if ( abs(cpress).lt. minc) cycle
c call TanhCorrel(span,span,thvec,sumth,sumth2,
c & vx,nvars,ixlo,ixhi,iylo,iyhi,i,j,1,direction,arho,brho,crho)
c if ( abs(brho) .lt. minjump*abs(arho) ) cycle
c if ( abs(crho).lt. minc) cycle
c if ( crho*cpress .lt. 0.0 ) then
slow = max(i-dcspan,1)
shigh = min(i+dcspan+1,nx+1)
do m=slow,shigh
dcflag(m,j,direction) = CLES_SWITCH_WENO
enddo
c endif
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,iylo,iyhi,i,j,5,direction,
$ apress,bpress,cpress)
if ( abs(bpress) .lt. minjump*abs(apress) ) cycle
if ( abs(cpress).lt. minc) cycle
c call TanhCorrel(span,span-1,thvecbis,sumthbis,sumth2bis,
c & vx,nvars,ixlo,ixhi,iylo,iyhi,i,j,1,direction,arho,brho,crho)
c if ( abs(brho) .lt. minjump*abs(arho)) cycle
c if ( abs(crho).lt. minc) cycle
c if ( crho*cpress .lt. 0.0 ) then
slow = max(i-dcspan,1)
shigh = min(i+dcspan,nx+1)
do m=slow,shigh
dcflag(m,j,direction) = CLES_SWITCH_WENO
enddo
c endif
endif
enddo
enddo
else
c This is esentially (nghost - span)*2.
nspany = (1-iylo-span)*2
if ( nspany .lt. 2 ) then
print *, 'Number of ghost cells is too low for'//
$ ' non-linear projection flagging'
stop
endif
do i=1, nx, 1
do l=1-nspany, 1+ny*2+nspany, 1
if ( mod(l,2) .eq. 0 ) then
c Center of cell
j = 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,iylo,iyhi,i,j,5,direction,
$ apress,bpress,cpress)
if ( abs(bpress) .lt. minjump*abs(apress) ) cycle
if ( abs(cpress).lt. minc) cycle
c call TanhCorrel(span,span,thvec,sumth,sumth2,
c & vx,nvars,ixlo,ixhi,iylo,iyhi,i,j,1,direction,arho,brho,crho)
c if ( abs(brho) .lt. minjump*abs(arho) ) cycle
c if ( abs(crho).lt. minc) cycle
c if ( crho*cpress .lt. 0.0 ) then
slow = max(j-dcspan,1)
shigh = min(j+dcspan+1,ny+1)
do m=slow,shigh
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
c endif
else
c Wall of cell
j = 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,iylo,iyhi,i,j,5,direction,
$ apress,bpress,cpress)
if ( abs(bpress) .lt. minjump*abs(apress) ) cycle
if ( abs(cpress).lt. minc) cycle
c call TanhCorrel(span,span-1,thvecbis,sumthbis,sumth2bis,
c & vx,nvars,ixlo,ixhi,iylo,iyhi,i,j,1,direction,arho,brho,crho)
c if ( abs(brho) .lt. minjump*abs(arho)) cycle
c if ( abs(crho).lt. minc) cycle
c if ( crho*cpress .lt. 0.0 ) then
slow = max(j-dcspan,1)
shigh = min(j+dcspan,ny+1)
do m=slow,shigh
dcflag(i,m,direction) = CLES_SWITCH_WENO
enddo
c endif
endif
enddo
enddo
endif
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
c Flag the boundaries
if ( direction .eq. 1 ) then
do j=1, ny
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,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
do j=1, ny
do i=1, nx
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
c--------------------------------------------------------------------------
return
end
c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
subroutine TanhCorrel(span,span1,thvec,sumth,sumth2,
& vx,nvars,ixlo,ixhi,iylo,iyhi,i,j,nfield,direction,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,iylo,iyhi,i,j,nfield,direction
double precision thvec(-span:span1)
double precision vx(nvars,ixlo:ixhi,iylo:iyhi)
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
if ( direction .eq. 1 ) then
do m=-span,span1, 1
sumy = sumy + vx(nfield,i+m,j)
sumyt = sumyt + vx(nfield,i+m,j)*thvec(m)
enddo
else
do m=-span,span1, 1
sumy = sumy + vx(nfield,i,j+m)
sumyt = sumyt + vx(nfield,i,j+m)*thvec(m)
enddo
endif
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
if ( direction .eq. 1 ) then
do m=-span,span1, 1
dp = vx(nfield,i+m,j)-a
sum = sum + dp*thvec(m)
sumy = sumy + dp*dp
enddo
else
do m=-span,span1, 1
dp = vx(nfield,i,j+m)-a
sum = sum + dp*thvec(m)
sumy = sumy + dp*dp
enddo
endif
if ( abs(sumy) .lt. 1.0d-12 ) then
c = 0.0d0
else
c = sum/sqrt(sumy*sumth2)
endif
return
end