cdl_hllc.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\cdl_hllc.f
00002 !
00072                         SUBROUTINE CDL_HLLC
00073 !                       ********************
00074 !
00075      &(NS,NPTFR,NBOR,LIMPRO,XNEBOR,YNEBOR,
00076      & W,CE,FLUENT,FLUSORT,FLBOR,EPS,WINF,
00077      & G,HBOR,UBOR,VBOR)
00078 !
00079 !
00080 !***********************************************************************
00081 ! TELEMAC 2D VERSION 6.2                                     01/15/2012
00082 !***********************************************************************
00083 !
00084 !
00085 !    UA(1,IS) = H,  UA(2,IS)=U  ,UA(3,IS)=V
00086 !
00087 !
00088 !
00089 !
00090 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00091 !|  NS            |-->|  TOTAL NUMBER OF NODES
00092 !|  NPTFR         |-->|  TOTAL NUMBER OF BOUNDARY NODES
00093 !|  NBOR          |-->|  GLOBAL NUMBERS OF BOUNDARY POINTS
00094 !|  LIMPRO        |-->|  TYPES OF BOUNDARY CONDITION
00095 !|  XNEBOR,YNEBOR |-->|  UNIT OUTWARD NORMAL COMPONENTS AT BOUNDARY POINTS
00096 !|  KDIR          |-->|  CONVENTION FOR DIRICHLET POINTS
00097 !|  KNEU          |-->|  CONVENTION FOR NEUMANN POINTS
00098 !|  W             |-->|  UA(1,IS) = H,  UA(2,IS)=U  ,UA(3,IS)=V
00099 !|  CE            |<->|  FLUX
00100 !|  FLUENT,FLUSORT|<--|  IN AND OUT MASS FLUX
00101 !|  FLBOR         |<--|  IN AND OUT WATER MASS FLUX
00102 !|  EPS           |-->|  TOLERANCE FOR WATER DEPTH DIVISION
00103 !|  HBOR,UBOR,VBOR|-->|  PRESCRIBED H, U AND V GIVEN BY BORD
00104 !|  WINF          |-->|  PRESCRIBED BOUNDARY CONDITIONS
00105 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00106 !
00107       USE BIEF
00108       USE DECLARATIONS_TELEMAC,ONLY: KDIR,KENT,KDDL,KNEU,KENTU
00109       USE DECLARATIONS_TELEMAC2D,ONLY:NUMLIQ,LIUBOR,ENTET
00110       USE INTERFACE_TELEMAC2D, EX_CDL_HLLC => CDL_HLLC
00111 !
00112       IMPLICIT NONE
00113       INTEGER LNG,LU
00114       COMMON/INFO/LNG,LU
00115 !
00116 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00117 !
00118       INTEGER, INTENT(IN)             :: NS,NPTFR
00119       INTEGER, INTENT(IN)             :: NBOR(NPTFR),LIMPRO(NPTFR,6)
00120       DOUBLE PRECISION, INTENT(IN)    :: XNEBOR(2*NPTFR),YNEBOR(2*NPTFR)
00121       DOUBLE PRECISION, INTENT(IN)    :: UBOR(NPTFR),VBOR(NPTFR)
00122       DOUBLE PRECISION, INTENT(IN)    :: HBOR(NPTFR)
00123       DOUBLE PRECISION, INTENT(IN)    :: W(3,NS),EPS,G
00124       DOUBLE PRECISION, INTENT(INOUT) :: WINF(3,NPTFR)
00125       DOUBLE PRECISION, INTENT(INOUT) :: CE(NS,3),FLUENT,FLUSORT
00126       TYPE(BIEF_OBJ) , INTENT(INOUT)  :: FLBOR
00127 !
00128 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00129 !
00130       INTEGER IS,K,IDRY
00131 !
00132       DOUBLE PRECISION :: VNX,VNY,XNN,YNN,VNL
00133       DOUBLE PRECISION :: UNN,VNN,LAMBDA1,LAMBDA2
00134       DOUBLE PRECISION :: FLX(4),H1,U10,U1,V1
00135       DOUBLE PRECISION :: H2,U2,V2,UGN,AL,C1
00136       DOUBLE PRECISION :: HG,UG,VG,CG,UGTEMP
00137       DOUBLE PRECISION :: INFLOW,OUTFLOW,REGIME,FOURG
00138       DOUBLE PRECISION,PARAMETER ::XI=0.0D0
00139 ! TO CORRECT WHEN CONSIDERING TRACERS
00140       DOUBLE PRECISION,PARAMETER ::PSI1=0.0D0,PSI2=0.0D0
00141       LOGICAL                    ::ROT,DEJA
00142 !
00143       DEJA =.FALSE.
00144       ROT = .TRUE.
00145       FOURG = 4.D0*G
00146 !
00147 ! LOOP OVER BOUNDARY NODES
00148       IF(NPTFR.GT.0)THEN    ! FOR PARALLEL CASES
00149       DO K=1,NPTFR
00150         IS=NBOR(K)
00151 !
00152 ! INITIALIZATION
00153         FLUENT  = 0.D0
00154         FLUSORT = 0.0D0
00155         INFLOW  = 0.0D0
00156         OUTFLOW = 0.0D0
00157         FLX(1)  = 0.0D0
00158         FLX(2)  = 0.0D0
00159         FLX(3)  = 0.0D0
00160         FLX(4)  = 0.0D0
00161 ! INDICATOR FOR DRY CELLS
00162         IDRY=0
00163 !   NORMALIZED NORMAL
00164         XNN=XNEBOR(K)
00165         YNN=YNEBOR(K)
00166 !   NON NORMALIZED NORMAL
00167         VNX=XNEBOR(K+NPTFR)
00168         VNY=YNEBOR(K+NPTFR)
00169 !
00170         VNL=SQRT(VNX**2+VNY**2)
00171 !
00172         H1 = W(1,IS)
00173         IF(H1.GT.EPS)THEN
00174           U1 = W(2,IS)/H1
00175           V1 = W(3,IS)/H1
00176         ELSE
00177           U1 = 0.0D0
00178           V1 = 0.0D0
00179           IDRY=IDRY+1
00180         ENDIF
00181 !**************************************************
00182 !         WALL BOUNDARY
00183 !**************************************************
00184 !===============================
00185 !    SLIPPING CONDITION
00186 !===============================
00187 !
00188         IF(LIMPRO(K,1).EQ.KNEU) THEN
00189 
00190 ! FIRST METHOD: STRONG IMPOSITION
00191 !********************************
00192 !    CE.n = 0  is done in cdlproj
00193 ! DEFINITION OF THE GHOST STATE Ue
00194           H2=H1
00195 !         ROTATION
00196           U10 = U1
00197           U1  = XNN*U10+YNN*V1
00198           V1  =-YNN*U10+XNN*V1
00199 ! PUT NORMAL COMPONENT = 0
00200           U1 =  0.0D0
00201           U2 =  U1
00202           V2 =  V1
00203 ! INVERSE ROTATION
00204           U10 = U1
00205           U1  = -YNN*V1
00206           V1  =  XNN*V1
00207 !
00208           U2  = -YNN*V2
00209           V2  =  XNN*V2
00210 ! SECOND METHOD: WEAK IMPOSITION
00211 !********************************
00212 !DEFINITION OF THE GHOST STATE Ue
00213 !           H2 = H1
00214 ! INNER PRODUCT 2V.n
00215 !           U10 = 2.D0*(U1*XNN + V1*YNN)
00216 ! WEAK IMPOSITION: PUT Ve = V1-2(V1.n)n
00217 !           U2 = U1 - U10*XNN
00218 !           V2 = V1 - U10*YNN
00219 !
00220           CALL FLUX_HLLC(XI,H1,H2,U1,U2,V1,V2,PSI1,PSI2,
00221      &                 XNN,YNN,ROT,FLX)
00222           GOTO 100
00223 !
00224 !**************************************************
00225 !         LIQUID BOUNDARIES
00226 !**************************************************
00227         ELSEIF((LIMPRO(K,1).EQ.KDIR).OR.(LIMPRO(K,1).EQ.KDDL))THEN
00228 !         PREPARE COMPUTATION OF RIEMANN INVARIANTS
00229           IF(H1.LT.EPS)THEN
00230             UNN = 0.D0
00231             VNN = 0.D0
00232           ELSE
00233             UNN =  XNN*U1 + YNN*V1
00234             VNN = -YNN*U1 + XNN*V1
00235           ENDIF
00236 !===============================
00237 !    IF H IS IMPOSED
00238 !===============================
00239 !
00240           IF(LIMPRO(K,1).EQ.KDIR) THEN
00241 !
00242             HG = HBOR(K) ! THIS IS HG (GHOST STATE)
00243             CG = SQRT(G*HG)
00244             C1 = SQRT(G*H1)
00245             LAMBDA1 = UNN + C1 ! WE USE REAL H (H1) TO ASSESS THE REGIME
00246             LAMBDA2 = UNN - C1
00247             REGIME  = LAMBDA1*LAMBDA2
00248 !
00249             IF(REGIME.LT.0.D0.OR.UNN.LE.0.D0) THEN ! SUBCRITICAL REGIME OR ENTRY
00250               IF(HG.LT.EPS)THEN
00251                 UG = 0.D0 !  UG (GHOST)
00252                 VG = 0.D0 !  VG (GHOST)
00253                 IDRY = IDRY + 1
00254               ELSE
00255                 IF(REGIME.LT.0.D0) THEN !SUBCRITICAL
00256                   UG = UNN +2.D0*(C1-CG)
00257                   VG = VNN
00258 !                 INVERSE ROTATION
00259                   UGTEMP = UG
00260                   UG = XNN*UGTEMP - YNN*VG
00261                   VG = YNN*UGTEMP + XNN*VG
00262                 ELSE                            ! SUPERCRITICAL
00263                   IF(LIUBOR%I(K).EQ.KENTU.OR.LIUBOR%I(K).EQ.KENT) THEN ! IMPOSED INFLOW
00264                                                                    ! OR IMPOSED VELOCITY
00265                     UG = UBOR(K)  ! FORCING H TO HG AND U TO UG (SUPERCRITICAL)
00266                     VG = VBOR(K)
00267                   ELSE ! DATA MISSING
00268                     UG = 0.D0 !THIS IS LAKE AT REST (HYPOTHESIS)
00269                     VG = 0.D0
00270                     IF(.NOT.DEJA.AND.ENTET)THEN
00271                       IF(LNG.EQ.1) WRITE(LU,30) NUMLIQ%I(K)
00272                       IF(LNG.EQ.2) WRITE(LU,31) NUMLIQ%I(K)
00273                       DEJA=.TRUE.
00274                     ENDIF
00275                   ENDIF
00276                 ENDIF
00277               ENDIF
00278             ELSE !THIS IS A SUPERCRITICAL OUTFLOW (NO NEED FOR GIVEN H)
00279               IF(.NOT.DEJA.AND.ENTET)THEN
00280                 IF(LNG.EQ.1) WRITE(LU,60) NUMLIQ%I(K)
00281                 IF(LNG.EQ.2) WRITE(LU,61) NUMLIQ%I(K)
00282                 DEJA=.TRUE.
00283 !               NO CONTRIBUTION
00284             ENDIF
00285           ENDIF
00286 !
00287           GOTO 90 ! TO COMPUTE THE FLUX
00288 !==================================
00289 !    IF GIVEN VELOCITY OR DISCHARGE
00290 !==================================
00291 !
00292         ELSEIF(LIMPRO(K,2).EQ.KDIR)THEN     !   (LIUBOR%I(K).EQ.KENTU)THEN
00293 !
00294           UG = UBOR(K) ! GIVEN BY BORD
00295           VG = VBOR(K) ! GIVEN BY BORD
00296 !
00297           UGN =  XNN*UG + YNN*VG ! TO RETRIEVE NORMAL COMPONENT OF UG
00298 !          VGN = -YNN*UG + XNN*VG ! AND SO ON
00299 !
00300           HG = (UNN + 2.D0*SQRT(G*H1)-UGN)**2/FOURG
00301 !
00302           IF(HG.LE.EPS)THEN
00303             IDRY = IDRY + 1
00304           ENDIF
00305 !
00306           GOTO 90 ! TO COMPUTE THE FLUX
00307 !=========================================================
00308 !    IF GIVEN DISCHARGE :CONSIDERED IN THE PREVIOUS CASES
00309 !    THANKS TO SUBROUTINE DEBIMP WHICH TRANSFORMS FLOWRATES
00310 !    TO H (AT TN, WHICH CAN BE IMPROVED) AND VELOCITY
00311 !==========================================================
00312 !         ELSE IF(LIUBOR%I(K).EQ.KENT)THEN
00313 !
00314           ENDIF
00315         ENDIF
00316 !
00317 90      CONTINUE
00318 !       COMPUTE THE FLUX
00319         IF(IDRY.LT.2)THEN
00320 !        AT LEAST ONE WET CELL
00321           CALL FLUX_HLLC(XI,H1,HG,U1,UG,V1,VG,PSI1,PSI2,
00322      &                    XNN,YNN,ROT,FLX)
00323         ENDIF
00324 !
00325 !       FINAL BALANCE
00326         OUTFLOW  = FLX(1)*VNL
00327         IF(UNN.LE.0D0)THEN
00328           FLUENT = FLUENT + OUTFLOW
00329         ELSE
00330           FLUSORT = FLUSORT + OUTFLOW
00331         ENDIF
00332         FLBOR%R(K) = OUTFLOW
00333 !
00334 100     CONTINUE
00335 !
00336         CE(IS,1)  = CE(IS,1) - VNL*FLX(1)
00337         CE(IS,2)  = CE(IS,2) - VNL*FLX(2)
00338         CE(IS,3)  = CE(IS,3) - VNL*FLX(3)
00339       ENDDO
00340 !
00341 !      PROJECTION ON THE NORMAL TO ELIMINATE THE TANGENT COMPONENT
00342 !      ONLY FOR LIQUID BOUNDARY
00343       DO K=1,NPTFR
00344         IS=NBOR(K)
00345         IF(LIMPRO(K,1).NE.KNEU) THEN
00346           !NORMALIZED NORMAL
00347           XNN=XNEBOR(K)
00348           YNN=YNEBOR(K)
00349           UGN =  XNN*CE(IS,2) + YNN*CE(IS,3) ! TO RETRIEVE NORMAL COMPONENET OF UG
00350 !         VGN =  0.D0  ! PUT TANGENTIAL COMPENENT =0
00351 !         INVERSE ROTATION
00352           CE(IS,2) = XNN*UGN
00353           CE(IS,3) = YNN*UGN
00354         ENDIF
00355       ENDDO
00356 !
00357       ENDIF ! PARALLEL CASES
00358 !
00359 30    FORMAT(1X,'CDL_HLLC: ATTENTION SUR LA FRONTIERE ',1I6,/,1X,
00360      & '          ENTREE TORRENTIELLE        ',/,1X,
00361      & '          ET DEBIT NON FOURNI',/,1X)
00362 31    FORMAT(1X,'CDL_HLLC: WARNING, LIQUID BOUNDARY ',1I6,/,1X,
00363      & '          SUPERCRITICAL INLET        ',/,1X,
00364      & '          AND NO DISCHARGE PROVIDED',/,1X)
00365 !
00366 60    FORMAT(1X,'CDL_HLLC: ATTENTION SUR LA FRONTIERE ',1I6,/,1X,
00367      & '          SORTIE TORRENTIELLE ET DONC       ',/,1X,
00368      & '          CONDITION AUX LIMITES NON IMPOSEE',/,1X)
00369 61    FORMAT(1X,'CDL_HLLC: WARNING, LIQUID BOUNDARY ',1I6,/,1X,
00370      & '          SUPERCRITICAL OUTLET        ',/,1X,
00371      & '          DESIRED BOUNDARY CONDITION MAY BE UNSATISFIED')
00372 !
00373 !-----------------------------------------------------------------------
00374 !
00375       RETURN
00376       END
00377 !**********************************************************************

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0