c xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx c c Curvature 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 double precision curvP, curvRho, cP, cRho integer span double precision CurvThressholdP, CurvThressholdRho CurvThressholdP = 0.0d0 CurvThressholdRho = 0.0d0 span = 3 ! ---- test criteria at each point if ( direction .eq. 1 ) then cP = CurvThressholdP * dx * dx * 0.25d0 cRho = CurvThressholdRho * dx * dx * 0.25d0 do k = 1, nz, 1 do j = 1, ny, 1 DO i = 1, nx+1, 1 curvP = (vx(5,i-1,j,k)-2.*vx(5,i,j,k)+vx(5,i+1,j,k)) curvP = curvP/(vx(5,i-1,j,k)+2.*vx(5,i,j,k) $ +vx(5,i+1,j,k)) curvrho = (vx(1,i-1,j,k)-2.*vx(1,i,j,k)+vx(1,i+1,j,k)) curvrho = curvrho/(vx(1,i-1,j,k)+2.*vx(1,i,j,k) $ +vx(1,i+1,j,k)) IF ((abs(curvP).gt. cP).and.(abs(curvrho).gt.cRho) $ .and.(curvP*curvrho.gt.0.0)) 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 cP = CurvThressholdP * dy * dy * 0.25d0 cRho = CurvThressholdRho * dy * dy * 0.25d0 do k = 1, nz, 1 do j = 1, ny+1, 1 do i = 1, nx, 1 curvP = (vx(5,i,j-1,k)-2.*vx(5,i,j,k)+vx(5,i,j+1,k)) curvP = curvP/(vx(5,i,j-1,k)+2.*vx(5,i,j,k) $ +vx(5,i,j+1,k)) curvrho = (vx(1,i,j-1,k)-2.*vx(1,i,j,k) $ +vx(1,i,j+1,k)) curvrho = curvrho/(vx(1,i,j-1,k)+2.*vx(1,i,j,k) $ +vx(1,i,j+1,k)) IF ((curvP.gt. cP).and.(curvrho.gt.cRho) $ .and.(curvP*curvrho.gt.0.0)) 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 cP = CurvThressholdP * dz * dz * 0.25d0 cRho = CurvThressholdRho * dz * dz * 0.25d0 do k = 1, nz+1, 1 do j = 1, ny, 1 DO i = 1, nx, 1 curvP = (vx(5,i,j,k-1)-2.*vx(5,i,j,k)+vx(5,i,j,k+1)) curvP = curvP/(vx(5,i,j,k-1)+2.*vx(5,i,j,k) $ +vx(5,i,j,k+1)) curvrho = (vx(1,i,j,k-1)-2.*vx(1,i,j,k) $ +vx(1,i,j,k+1)) curvrho = curvrho/(vx(1,i,j,k-1)+2.*vx(1,i,j,k) $ +vx(1,i,j,k+1)) IF ((curvP.gt. cP).and.(curvrho.gt.cRho) $ .and.(curvP*curvrho.gt.0.0)) 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