flux_hllc.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac2d\flux_hllc.f
00002 !
00043                          SUBROUTINE FLUX_HLLC
00044 !                       ***********************
00045      &(XI,H1,H2,U1,U2,V1,V2,PSI1,PSI2,
00046      & XNN,YNN,ROT,HLLCFLX)
00047 !
00048 !***********************************************************************
00049 ! TELEMAC 2D VERSION 7.0                                         R. ATA
00050 !
00051 !***********************************************************************
00052 !BRIEF
00053 !
00054 !     FUNCTION  : SUBROUTINE COMPUTES HLLC FLUX: THREE HYDRODYNAMICAL
00055 !                 COMPENENTS + TRACER TRANSPORT
00056 !      SEE TORO: SHOCK CAPTURING METHODS FOR FREE
00057 !            SURFACE FLOWS (WILEY 2005)
00058 !
00059 !HISTORY  RIADH ATA (EDF R&D-LNHE)
00060 !+        07/15/2012
00061 !+        V6P2
00062 !+
00063 !
00064 !HISTORY  RIADH ATA (EDF R&D-LNHE)
00065 !+        03/20/2013
00066 !+        V6P3
00067 !+  OPTIMIZATION OF THE CODE
00068 !+  AVOID DIVISION BY 0
00069 !
00070 !HISTORY  RIADH ATA (EDF R&D-LNHE)
00071 !+        10/6/2013
00072 !+        V6P3
00073 !+  BUG FIXED IN COMPUTING U*
00074 !+  THANKS TO L. STADLER (BAW)
00075 !
00076 !
00077 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00078 ! |  FLX           | <--|  FLUX COMPONENTS AT THE INTERFACE            |
00079 ! |  H1,H2         | -->|  LEFT AND RIGHT WATER DEPTHS                 |
00080 ! |  PSI1,PSI2     | -->|  LEFT AND RIGHT TRACER DENSITIES             |
00081 ! |  ROT           | -->|  EXECUTE FINAL ROTATION OR NO                |
00082 ! |  U1,U2         | -->|  LEFT AND RIGHT VELOCITY X-COMPONENTS        |
00083 ! |  V1,V2         | -->|  LEFT AND RIGHT VELOCITY Y-COMPONENTS        |
00084 ! |  XNN,YNN       | -->|  X AND Y COMPONENT OF THE OUTWARD NORMAL     |
00085 ! ______________________________________________________________________
00086 !
00087 !  MODE: -->(UNCHANGEABLE INPUT),<--(OUTPUT),<-->(CHANGEABLE INPUT)
00088 !-----------------------------------------------------------------------
00089 !  CALLING SUBROUTINE FLUX_WAF OR FLUX_HLLC OR FLUXZZ
00090 !
00091 !***********************************************************************
00092 !
00093       USE BIEF
00094 !
00095       IMPLICIT NONE
00096       INTEGER LNG,LU
00097       COMMON/INFO/LNG,LU
00098 !
00099 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00100 !
00101       DOUBLE PRECISION, INTENT(IN)    :: H1,H2,U1,U2,PSI1,PSI2
00102       DOUBLE PRECISION, INTENT(IN)    :: XI,V1,V2,XNN,YNN
00103       LOGICAL, INTENT(IN)             :: ROT
00104       DOUBLE PRECISION, INTENT(INOUT) :: HLLCFLX(4)
00105 !***********************************************************************
00106 !
00107       INTEGER                         :: I,SPY
00108       DOUBLE PRECISION, PARAMETER     :: G = 9.81D0
00109       DOUBLE PRECISION                :: HL,UL,VL,PSI_L
00110       DOUBLE PRECISION                :: HR,UR,VR,PSI_R
00111       DOUBLE PRECISION                :: AL,AR,HSTAR,USTAR
00112       DOUBLE PRECISION                :: PQL,PQR,SL,SR
00113       DOUBLE PRECISION                :: QSTARL(4),QSTARR(4)
00114       DOUBLE PRECISION                :: QL(4),QR(4),FL(4),FR(4)
00115       DOUBLE PRECISION                :: FSTARL(4),FSTARR(4)
00116 !
00117       DOUBLE PRECISION                :: GSUR2,EPS,DENOM
00118       DOUBLE PRECISION                :: FLU2X,FLU2Y
00119       DOUBLE PRECISION                :: U0,POND,SSTAR
00120       DOUBLE PRECISION                :: FLX(4)
00121 !-----------------------------------------------------------------------
00122       EPS   = 1.E-6
00123       GSUR2 = G/2.0D0
00124       SPY   = 0
00125       PQL   = 0.0D0
00126       PQR   = 0.0D0
00127       USTAR = 0.0D0
00128       HSTAR = 0.0D0
00129       AL    = 0.0D0
00130       AR    = 0.0D0
00131 !***********************************************************************
00132 ! INITIALIZATION OF FLX AND HLLCFLX
00133       DO I=1,4
00134         FLX(I)     = 0.0D0
00135         HLLCFLX(I) = 0.0D0
00136       ENDDO
00137 !
00138 !-----------------------------------------------------------------------
00139 ! DEPTHS, VELOCITIES, TRACERS
00140       HL    = H1
00141       UL    = U1
00142       VL    = V1
00143       PSI_L = PSI1
00144 !
00145       HR    = H2
00146       UR    = U2
00147       VR    = V2
00148       PSI_R = PSI2
00149 !
00150 ! ROTATION
00151 !
00152       U0 = UL
00153       UL  = XNN*U0+YNN*VL
00154       VL  =-YNN*U0+XNN*VL
00155 !
00156       U0 = UR
00157       UR  = XNN*U0+YNN*VR
00158       VR  =-YNN*U0+XNN*VR
00159 !
00160 ! CASE WITH DRY LEFT AND RIGHT
00161       IF(HL.LT.EPS.AND.HR.LT.EPS)GOTO 20
00162 !
00163 ! CELERITIES
00164 !
00165       AL = SQRT(G*HL)
00166       AR = SQRT(G*HR)
00167 ! STAR VARIABLES
00168       HSTAR = 0.5D0*(HL+HR)-0.25D0*(UR-UL)*(HL+HR)/(AL+AR)
00169 !RA BUG FIXED WHEN COMPUTING U STAR
00170 !       USTAR = 0.5D0*(UL+UR)-0.25D0*(HR-HL)*(AL+AR)/(HL+HR)
00171       USTAR = 0.5D0*(UL+UR)-       (HR-HL)*(AL+AR)/(HL+HR)
00172 ! COMPUTE PQL AND PQR:
00173 ! IT WILL DEPEND IF WE ARE IN PRESENCE OF SHOCK OR RAREFACTION WAVE
00174       IF(HSTAR.LT.HL)THEN
00175 !       RAREFACTION
00176         PQL = 1.0D0
00177       ELSE
00178 !       SHOCK
00179         IF(HL.GT.EPS)THEN
00180           PQL = SQRT(0.5D0*(HSTAR+HL)*HSTAR/HL**2)
00181         ELSE
00182           PQL = 0.0D0
00183         ENDIF
00184       ENDIF
00185       IF(HSTAR.LT.HR)THEN
00186 !       RAREFACTION
00187         PQR = 1.0D0
00188       ELSE
00189 !       SHOCK
00190         IF(HR.GT.EPS)THEN
00191           PQR = SQRT(0.5D0*(HSTAR+HR)*HSTAR/HR**2)
00192         ELSE
00193           PQR = 0.0D0
00194         ENDIF
00195       ENDIF
00196 !
00197 20    CONTINUE
00198 !
00199 ! COMPUTE SL, SR AND SSTAR  (WE CONSIDER DRY CASES)
00200       IF(HL.GT.EPS)THEN
00201         SL = UL-AL*PQL
00202       ELSE
00203         SL = UR - 2.0D0*AR
00204         SR = UR + AR
00205 ! RA+SP: USE OF ANALYTICAL FORMULA FOR SSTAR
00206 !       SSTAR = SL
00207         GOTO 35
00208 ! RA+SP: USE OF ANALYTICAL FORMULA FOR SSTAR
00209 !        SSTAR = SL
00210         GOTO 35
00211       ENDIF
00212 !
00213       IF(HR.GT.EPS)THEN
00214         SR = UR + AR*PQR
00215       ELSE
00216         SL = UL - AL
00217         SR = UL + 2.0D0*AL
00218 ! RA+SP: USE OF ANALYTICAL FORMULA FOR SSTAR
00219 !       SSTAR = SR
00220         GOTO 35
00221 ! RA+SP: USE OF ANALYTICAL FORMULA FOR SSTAR
00222 !         SSTAR = SR
00223         GOTO 35
00224       ENDIF
00225 !RA      SSTAR = USTAR
00226 35    CONTINUE
00227       DENOM = HR*(UR-SR)-HL*(UL-SL)
00228       IF(ABS(DENOM).LT.EPS)THEN
00229         SSTAR = USTAR
00230       ELSE
00231         SSTAR = (SL*HR*(UR-SR)-SR*HL*(UL-SL))/DENOM
00232       ENDIF
00233 !END RA
00234 !
00235 ! COMPUTE QL AND QR
00236       QL(1)     = HL
00237       QL(2)     = HL*UL
00238       QL(3)     = HL*VL
00239       QL(4)     = HL*PSI_L
00240 !
00241       QR(1)     = HR
00242       QR(2)     = HR*UR
00243       QR(3)     = HR*VR
00244       QR(4)     = HR*PSI_R
00245 
00246 ! COMPUTE QSTARL AND QSTARR
00247       IF(ABS(SL-SSTAR).GT.EPS)THEN
00248         POND = HL*(SL-UL)/(SL-SSTAR)
00249       ELSE
00250         POND = 0.0D0
00251       ENDIF
00252       QSTARL(1) = POND
00253       QSTARL(2) = POND*SSTAR
00254       QSTARL(3) = POND*VL
00255       QSTARL(4) = POND*PSI_L
00256 !
00257       IF(ABS(SR-SSTAR).GT.EPS)THEN
00258         POND = HR*(SR-UR)/(SR-SSTAR)
00259       ELSE
00260         POND = 0.0D0
00261       ENDIF
00262       QSTARR(1) = POND
00263       QSTARR(2) = POND*SSTAR
00264       QSTARR(3) = POND*VR
00265       QSTARR(4) = POND*PSI_R
00266 !
00267 ! COMPUTE FL AND FR
00268 !
00269       FL(1)     = HL*UL
00270       FL(2)     = HL*UL**2 +GSUR2*HL**2
00271       FL(3)     = HL*UL*VL
00272       FL(4)     = HL*UL*PSI_L
00273 !
00274       FR(1)     = HR*UR
00275       FR(2)     = HR*UR**2 +GSUR2*HR**2
00276       FR(3)     = HR*UR*VR
00277       FR(4)     = HR*UR*PSI_R
00278 !
00279 ! COMPUTE FSTARL SFTARR
00280       FSTARL(1) = FL(1) + SL*(QSTARL(1)-QL(1))
00281       FSTARL(2) = FL(2) + SL*(QSTARL(2)-QL(2))
00282       FSTARL(3) = FL(3) + SL*(QSTARL(3)-QL(3))
00283       FSTARL(4) = FL(4) + SL*(QSTARL(4)-QL(4))
00284 !
00285       FSTARR(1) = FR(1) + SR*(QSTARR(1)-QR(1))
00286       FSTARR(2) = FR(2) + SR*(QSTARR(2)-QR(2))
00287       FSTARR(3) = FR(3) + SR*(QSTARR(3)-QR(3))
00288       FSTARR(4) = FR(4) + SR*(QSTARR(4)-QR(4))
00289 ! AND FINALLY THE HLLC FLUX (BEFORE ROTATION)
00290       IF(XI.LT.SL)THEN
00291         FLX(1) = FL(1)
00292         FLX(2) = FL(2)
00293         FLX(3) = FL(3)
00294         FLX(4) = FL(4)
00295         SPY = 1
00296       ELSEIF(XI.LT.SSTAR.AND.XI.GT.SL) THEN
00297         FLX(1) = FSTARL(1)
00298         FLX(2) = FSTARL(2)
00299         FLX(3) = FSTARL(3)
00300         FLX(4) = FSTARL(4)
00301         SPY = 1
00302       ELSEIF(XI.GT.SSTAR.AND.XI.LT.SR) THEN
00303         FLX(1) = FSTARR(1)
00304         FLX(2) = FSTARR(2)
00305         FLX(3) = FSTARR(3)
00306         FLX(4) = FSTARR(4)
00307         SPY = 1
00308       ELSE
00309         FLX(1) = FR(1)
00310         FLX(2) = FR(2)
00311         FLX(3) = FR(3)
00312         FLX(4) = FR(4)
00313         SPY = 1
00314       ENDIF
00315       IF(SPY.EQ.0)THEN
00316         WRITE(LU,*)'ERROR IN HLLC FLUX ESTIMATION (FLUX_HLLC.F)'
00317         CALL PLANTE(1)
00318         STOP
00319       ENDIF
00320 !
00321 ! INVERSE ROTATION AND FINAL FLUX
00322 !
00323       IF(ROT)THEN
00324         FLU2X  = XNN*FLX(2) - YNN*FLX(3)
00325         FLU2Y  = YNN*FLX(2) + XNN*FLX(3)
00326 !
00327         HLLCFLX(1) = FLX(1)
00328         HLLCFLX(2) = FLU2X
00329         HLLCFLX(3) = FLU2Y
00330         HLLCFLX(4) = FLX(4)
00331       ELSE
00332 ! IN THIS CASE, NO ROTATION
00333 !
00334 ! FINAL FLUX
00335 !
00336         HLLCFLX(1) = FLX(1)
00337         HLLCFLX(2) = FLX(2)
00338         HLLCFLX(3) = FLX(3)
00339         HLLCFLX(4) = FLX(4)
00340       ENDIF
00341 !
00342 !-----------------------------------------------------------------------
00343 !
00344       RETURN
00345       END

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