walldist.f

Go to the documentation of this file.

!                   *******************
                    SUBROUTINE walldist
!                   *******************
     &(w_dist,t11,t12,t13,t14,t15,flbor,tb,am1,am2,s,
     & liubor,ielmnu,nptfr,mesh)

!***********************************************************************
! TELEMAC2D   V7P3                                             11/2016
!***********************************************************************
!
!brief    COMPUTES THE DISTANCE TO THE CLOSEST WALL FOR
!              SPALART ALLMARAS MODEL-2D COMPUTATION.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AM1,AM2        |<--| WORKING MATRICES
!| LIUBOR         |-->| BOUNDARY CONDITION ON VELOCITIES
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| S              |<--| VOID STRUCTURE
!| T12..FLBOR     |<--| WORKING ARRAYS
!| WDIST          |<--| DISTANCE FROM THE CLOSEST WALL
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!***********************************************************************
      USE bief, ex_walldist => walldist
      USE declarations_special
      USE declarations_parallel
      USE interface_parallel
      USE declarations_telemac, ONLY: kadh,klog,kneu,kdir
!
      IMPLICIT NONE
      INTRINSIC sqrt
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER        , INTENT(IN   ) :: IELMNU,NPTFR
      TYPE(bief_obj) , INTENT(IN   ) :: LIUBOR
      TYPE(bief_obj) , INTENT(INOUT) :: W_DIST
      TYPE(bief_obj) , INTENT(INOUT) :: T11,T12,T13,T14,T15,FLBOR,TB
      TYPE(bief_obj) , INTENT(INOUT) :: AM1,AM2,S
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      LOGICAL :: SOLVER_INFO
      INTEGER :: IPTFR,IP,IELWD,IERR,DEBUG
      DOUBLE PRECISION  :: GRAD_D
      INTEGER, ALLOCATABLE :: LIMDIST(:)
      TYPE(slvcfg) :: SLVDIST
!
!-----------------------------------------------------------------------
      ALLOCATE(limdist(nptfr))
      ierr = 0
      debug= 0
!     SETUP LINEAR SOLVER CONFIGURATION STRUCTURE
      slvdist%SLV = 1
      slvdist%PRECON = 2
      slvdist%KRYLOV = 0
      slvdist%NITMAX = 1000
      slvdist%EPS = 1.d-8
      IF(debug.GT.0) THEN
        solver_info = .true.
      ELSE
        solver_info = .false.
      ENDIF
!
      ielwd = ielmnu ! ELEMENTS USED FOR S-A
            WRITE(lu,*) 'TELEMAC2D: 31 '
!
      DO ip=1,w_dist%DIM1
        t11%R(ip) = 1.d0 ! VISC
        t12%R(ip) = 1.d0 ! RHS SOURCE TERM
        t13%R(ip) = 0.d0 ! WALL DIST
        t14%R(ip) = 0.d0 ! WORK
        t15%R(ip) = 0.d0 ! BCS
      ENDDO
!
      DO iptfr=1,nptfr
!
        IF(liubor%I(iptfr).EQ.klog.OR.liubor%I(iptfr).EQ.kadh) THEN
          limdist(iptfr) = kdir
          flbor%R(iptfr)   = 0.d0
        ELSE
          limdist(iptfr) = kneu
        ENDIF
!
      ENDDO
!
!     CREATES THE MASS MATRIX (AM1)
!
      CALL matrix(am1,'M=N     ','MATMAS          ',ielwd,ielwd,
     &            1.d0,s,s,s,s,s,s,mesh,.false.,s)
!
!     CREATES THE DIFFUSION MATRIX (AM2)
!
      CALL matrix(am2,'M=N     ','MATDIF          ',ielwd,ielwd,
     &            1.d0,s,s,s,t11,t11,t11,mesh,.false.,s)
!
!     PREPARES RHS
      CALL matvec('X=AY    ',t13,am1,t12,1.d0,mesh)
!
!     IMPOSES DIRICHLET CONDITIONS
!
!     BE AWARE: T1 AND T2 ARE ERASED BY DIRICH VIA TB
      IF(nptfr.GT.0) THEN
        CALL dirich(t14,am2,t13,flbor,limdist,
     &              tb,mesh,kdir,.false.,s)
      ENDIF
!
!     SOLVES THE LINEAR SYSTEM
!
!     BE AWARE: T1 -> T7 ARE ERASED BY SOLVE VIA TB
      IF (debug.GT.0) WRITE(lu,*) 'SOLVING WALL DISTANCE'
      CALL solve(t14,am2,t13,tb,slvdist,solver_info,mesh,am2)
!
!     COMPUTES VERTEX GRADIENT
      CALL vector(t11,'=','GRADF          X',ielwd,
     &            1.d0,t14,s,s,s,s,s,mesh,.false.,s,asspar=.true.)
      CALL vector(t12,'=','GRADF          Y',ielwd,
     &            1.d0,t14,s,s,s,s,s,mesh,.false.,s,asspar=.true.)
      CALL vector(t15 , '=' , 'MASBAS          ' , ielwd ,
     &            1.d0,s,s,s,s,s,s,
     &            mesh,.false.,s,asspar=.true.)
!     NORMALIZES GRADIENT
      CALL os('X=Y/Z   ', x=t11, y=t11, z=t15)
      CALL os('X=Y/Z   ', x=t12, y=t12, z=t15)
!
!     CALCULATES THE WALL DISTANCE
!
      DO ip=1,w_dist%DIM1
        grad_d = t11%R(ip)**2.d0 + t12%R(ip)**2.d0
        IF(grad_d + 2.d0*t14%R(ip).GT.0.d0.AND.grad_d.GT.0) THEN
          w_dist%R(ip) = -sqrt(grad_d) + sqrt(grad_d + 2.d0*t14%R(ip))
        ELSE
          w_dist%R(ip) = 0.d0
          ierr = ierr + 1
        ENDIF
      ENDDO
      IF(ncsize.GT.0)THEN
        ierr = p_sum(ierr)
      ENDIF
!
      IF(ierr.NE.0)THEN
        WRITE(lu,*) 'WALLDIST (BIEF): WARNING WITH THE WALL DISTANCE '
        WRITE(lu,*) '           COMPUTATION. NEGATIVE DETERMINANT '
        WRITE(lu,*) '           ON ', ierr, 'NODES'
        WRITE(lu,*) '           (PROBABLY OVER CONSTRAINT NODES)'
      ENDIF
      DEALLOCATE(limdist)
!
!-----------------------------------------------------------------------
!
      END SUBROUTINE