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 double precision sum, sumth, sumth2, sumy, sumyt, th double precision x, d, det, a, b, c, dp integer span, n, imin, imax parameter (span=3) double precision thvec(-span:span) c double precision minc, minjump, dwidth c thickness of the profile in terms of grid sizes c a value of 1 means a discontinuity of 6 grids points c because of the thickness of the tanh() c dwidth = 0.75d0 c minimum value of a correlation coefficient considered c to be a good match c minc = 0.9d0 c to discriminate small oscillations in the density that c have no meaning, we stipulate that any jump in density c that is below 2% of the mean value is not a jump. c minjump = 0.2d-1 c ---- test criteria at each point sumth = 0.0d0 sumth2 = 0.0d0 do m=-span,span, 1 x = m/dwidth th = tanh(x) thvec(m) = th sumth = sumth + th sumth2 = sumth2 + th*th enddo n = span*2+1 imin = max(ixlo+span,1) imax = min(ixhi-span,nx) DO i = imin, imax, 1 sumy = 0.0d0 sumyt = 0.0d0 do m=-span,span, 1 sumy = sumy + vx(1,i+m) sumyt = sumyt + vx(1,i+m)*thvec(m) enddo c Optimal values det = n*sumth2-sumth*sumth a = (sumy*sumth2-sumth*sumyt)/det b = (n*sumyt-sumth*sumy)/det c some tests conditions if ( a .lt. 0.0d0 ) cycle if ( a + b .lt. 0.0d0 ) cycle if ( a - b .lt. 0.0d0 ) cycle if ( abs(b) .lt. minjump*a ) cycle c Compute correlation sum = 0.0d0 sumy = 0.0d0 do m=-span,span, 1 dp = vx(1,i+m)-a sum = sum + dp*thvec(m) sumy = sumy + dp*dp enddo if ( abs(sumy) .lt. 1.0d-12 ) cycle c = sum/sqrt(sumy*sumth2) if ( abs(c) .gt. minc ) then slow = max(i-span,1) shigh = min(i+1+span,nx+1) do m=slow,shigh dcflag(m,direction) = CLES_SWITCH_WENO enddo endif enddo return end