flux_ef_vf.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\flux_ef_vf.f
00002 !
00076                      SUBROUTINE FLUX_EF_VF
00077 !                    *********************
00078 !
00079      &(FLOW,PHIEL,NSEG,NELEM,ELTSEG,ORISEG,IKLE,INIFLO,IOPT,FN)
00080 !
00081 !***********************************************************************
00082 ! BIEF   V7P0                                   21/08/2010
00083 !***********************************************************************
00084 !
00085 !
00086 !
00087 !
00088 !
00089 !
00090 !
00091 !
00092 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00093 !| ELTSEG         |-->| SEGMENTS OF EVERY TRIANGLE
00094 !| FLOW           |<--| FLUXES PER SEGMENTS
00095 !| FN             |-->| OPTIONAL ARGUMENT FOR PSI SCHEME
00096 !| IKLE           |-->| CONNECTIVITY TABLE
00097 !| INIFLO         |-->| IF(YES) FLOW WILL BE INITIALISED AT 0
00098 !| IOPT           |-->| OPTION FOR THE CONSTANT PER ELEMENT
00099 !| NELEM          |-->| NUMBER OF ELEMENTS
00100 !| NSEG           |-->| NUMBER OF SEGMENTS
00101 !| ORISEG         |-->| ORIENTATION OF SEGMENTS OF EVERY TRIANGLE
00102 !| PHIEL          |-->| PER ELEMENT, FLUXES LEAVING POINTS
00103 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00104 !
00105       USE BIEF, EX_FLUX_EF_VF => FLUX_EF_VF
00106 !
00107       IMPLICIT NONE
00108       INTEGER LNG,LU
00109       COMMON/INFO/LNG,LU
00110 !
00111 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00112 !
00113       INTEGER, INTENT(IN)                  :: NSEG,IOPT,NELEM
00114       INTEGER, INTENT(IN)                  :: ELTSEG(NELEM,3)
00115       INTEGER, INTENT(IN)                  :: ORISEG(NELEM,3)
00116       INTEGER, INTENT(IN)                  :: IKLE(NELEM,3)
00117       DOUBLE PRECISION, INTENT(INOUT)      :: FLOW(NSEG)
00118       DOUBLE PRECISION, INTENT(IN)         :: PHIEL(NELEM,3)
00119       LOGICAL, INTENT(IN)                  :: INIFLO
00120       TYPE(BIEF_OBJ), INTENT(IN), OPTIONAL :: FN
00121 !
00122 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00123 !
00124       INTEGER IELEM,ISEG,I1,I2,I3
00125       DOUBLE PRECISION A1,A2,A3,K1,K2,K3,THIRD,CSTE,F12,F23,F31,F1,F2,F3
00126       DOUBLE PRECISION FN1,FN2,FN3,F21,F32,F13
00127       DOUBLE PRECISION BETA1FI,BETA2FI,BETA3FI,FI
00128 !
00129       INTRINSIC ABS,MIN,MAX
00130 !
00131       THIRD=1.D0/3.D0
00132 !
00133 !-----------------------------------------------------------------------
00134 !
00135 !     INITIALISES FLOW TO 0.D0
00136 !
00137       IF(INIFLO) THEN
00138         DO ISEG = 1,NSEG
00139           FLOW(ISEG) = 0.D0
00140         ENDDO
00141       ENDIF
00142 !
00143 !-----------------------------------------------------------------------
00144 !
00145       IF(IOPT.EQ.-1) THEN
00146 !
00147 !-----------------------------------------------------------------------
00148 !
00149 !     FLUXES ALREADY COMPUTED BEFORE CALLING THIS SUBROUTINE
00150 !     THEY ARE JUST ASSEMBLED HERE
00151 !
00152       DO IELEM = 1,NELEM
00153 !       SEGMENT 1
00154         ISEG  = ELTSEG(IELEM,1)
00155         IF(ORISEG(IELEM,1).EQ.1) THEN
00156           FLOW(ISEG) = FLOW(ISEG) + PHIEL(IELEM,1)
00157         ELSE
00158           FLOW(ISEG) = FLOW(ISEG) - PHIEL(IELEM,1)
00159         ENDIF
00160 !       SEGMENT 2
00161         ISEG  = ELTSEG(IELEM,2)
00162         IF(ORISEG(IELEM,2).EQ.1) THEN
00163           FLOW(ISEG) = FLOW(ISEG) + PHIEL(IELEM,2)
00164         ELSE
00165           FLOW(ISEG) = FLOW(ISEG) - PHIEL(IELEM,2)
00166         ENDIF
00167 !       SEGMENT 3
00168         ISEG  = ELTSEG(IELEM,3)
00169         IF(ORISEG(IELEM,3).EQ.1) THEN
00170           FLOW(ISEG) = FLOW(ISEG) + PHIEL(IELEM,3)
00171         ELSE
00172           FLOW(ISEG) = FLOW(ISEG) - PHIEL(IELEM,3)
00173         ENDIF
00174       ENDDO
00175 !
00176 !-----------------------------------------------------------------------
00177 !
00178       ELSEIF(IOPT.EQ.0) THEN
00179 !
00180 !-----------------------------------------------------------------------
00181 !
00182 !     WITH NO CONSTANT
00183 !
00184       DO IELEM = 1,NELEM
00185         F1 = PHIEL(IELEM,1)
00186         F2 = PHIEL(IELEM,2)
00187         F3 = PHIEL(IELEM,3)
00188 !       SEGMENT 1
00189         ISEG  = ELTSEG(IELEM,1)
00190         IF(ORISEG(IELEM,1).EQ.1) THEN
00191           FLOW(ISEG) = FLOW(ISEG) + THIRD*(F1-F2)
00192         ELSE
00193           FLOW(ISEG) = FLOW(ISEG) - THIRD*(F1-F2)
00194         ENDIF
00195 !       SEGMENT 2
00196         ISEG  = ELTSEG(IELEM,2)
00197         IF(ORISEG(IELEM,2).EQ.1) THEN
00198           FLOW(ISEG) = FLOW(ISEG) + THIRD*(F2-F3)
00199         ELSE
00200           FLOW(ISEG) = FLOW(ISEG) - THIRD*(F2-F3)
00201         ENDIF
00202 !       SEGMENT 3
00203         ISEG  = ELTSEG(IELEM,3)
00204         IF(ORISEG(IELEM,3).EQ.1) THEN
00205           FLOW(ISEG) = FLOW(ISEG) + THIRD*(F3-F1)
00206         ELSE
00207           FLOW(ISEG) = FLOW(ISEG) - THIRD*(F3-F1)
00208         ENDIF
00209       ENDDO
00210 !
00211 !-----------------------------------------------------------------------
00212 !
00213       ELSEIF(IOPT.EQ.1) THEN
00214 !
00215 !-----------------------------------------------------------------------
00216 !
00217 !     MINIMISES MAX ( ABS(FLOW) )
00218 !
00219       DO IELEM = 1,NELEM
00220         F1 = PHIEL(IELEM,1)
00221         F2 = PHIEL(IELEM,2)
00222         F3 = PHIEL(IELEM,3)
00223         CSTE=-0.5D0*(MIN(F1-F2,F2-F3,F3-F1)+MAX(F1-F2,F2-F3,F3-F1))
00224 !       SEGMENT 1
00225         ISEG  = ELTSEG(IELEM,1)
00226         IF(ORISEG(IELEM,1).EQ.1) THEN
00227           FLOW(ISEG) = FLOW(ISEG) + THIRD*(F1-F2+CSTE)
00228         ELSE
00229           FLOW(ISEG) = FLOW(ISEG) - THIRD*(F1-F2+CSTE)
00230         ENDIF
00231 !       SEGMENT 2
00232         ISEG  = ELTSEG(IELEM,2)
00233         IF(ORISEG(IELEM,2).EQ.1) THEN
00234           FLOW(ISEG) = FLOW(ISEG) + THIRD*(F2-F3+CSTE)
00235         ELSE
00236           FLOW(ISEG) = FLOW(ISEG) - THIRD*(F2-F3+CSTE)
00237         ENDIF
00238 !       SEGMENT 3
00239         ISEG  = ELTSEG(IELEM,3)
00240         IF(ORISEG(IELEM,3).EQ.1) THEN
00241           FLOW(ISEG) = FLOW(ISEG) + THIRD*(F3-F1+CSTE)
00242         ELSE
00243           FLOW(ISEG) = FLOW(ISEG) - THIRD*(F3-F1+CSTE)
00244         ENDIF
00245       ENDDO
00246 !
00247 !-----------------------------------------------------------------------
00248 !
00249       ELSEIF(IOPT.EQ.2) THEN
00250 !
00251 !-----------------------------------------------------------------------
00252 !
00253 !     LEO POSTMA'S METHOD (EQUIVALENT TO FLUXES GIVEN BY N-SCHEME)
00254 !
00255       DO IELEM = 1,NELEM
00256         F1 = PHIEL(IELEM,1)
00257         F2 = PHIEL(IELEM,2)
00258         F3 = PHIEL(IELEM,3)
00259         A1 = ABS(F1)
00260         A2 = ABS(F2)
00261         A3 = ABS(F3)
00262         IF(A1.GE.A2.AND.A1.GE.A3) THEN
00263 !         ALL FLOW TO AND FROM NODE 1
00264           ISEG  = ELTSEG(IELEM,1)
00265           IF(ORISEG(IELEM,1).EQ.1) THEN
00266             FLOW(ISEG) = FLOW(ISEG) - F2
00267           ELSE
00268             FLOW(ISEG) = FLOW(ISEG) + F2
00269           ENDIF
00270           ISEG = ELTSEG(IELEM,3)
00271           IF(ORISEG(IELEM,3).EQ.2) THEN
00272             FLOW(ISEG) = FLOW(ISEG) - F3
00273           ELSE
00274             FLOW(ISEG) = FLOW(ISEG) + F3
00275           ENDIF
00276         ELSEIF(A2.GE.A1.AND.A2.GE.A3) THEN
00277 !         ALL FLOW TO AND FROM NODE 2
00278           ISEG = ELTSEG(IELEM,1)
00279           IF(ORISEG(IELEM,1).EQ.2) THEN
00280             FLOW(ISEG) = FLOW(ISEG) - F1
00281           ELSE
00282             FLOW(ISEG) = FLOW(ISEG) + F1
00283           ENDIF
00284           ISEG = ELTSEG(IELEM,2)
00285           IF(ORISEG(IELEM,2).EQ.1) THEN
00286             FLOW(ISEG) = FLOW(ISEG) - F3
00287           ELSE
00288             FLOW(ISEG) = FLOW(ISEG) + F3
00289           ENDIF
00290         ELSE
00291 !         ALL FLOW TO AND FROM NODE 3
00292           ISEG = ELTSEG(IELEM,2)
00293           IF(ORISEG(IELEM,2).EQ.2) THEN
00294             FLOW(ISEG) = FLOW(ISEG) - F2
00295           ELSE
00296             FLOW(ISEG) = FLOW(ISEG) + F2
00297           ENDIF
00298           ISEG = ELTSEG(IELEM,3)
00299           IF(ORISEG(IELEM,3).EQ.1) THEN
00300             FLOW(ISEG) = FLOW(ISEG) - F1
00301           ELSE
00302             FLOW(ISEG) = FLOW(ISEG) + F1
00303           ENDIF
00304         ENDIF
00305       ENDDO
00306 !
00307 !-----------------------------------------------------------------------
00308 !
00309       ELSEIF(IOPT.EQ.3.AND.PRESENT(FN)) THEN
00310 !
00311 !-----------------------------------------------------------------------
00312 !
00313 !     PSI-SCHEME
00314 !
00315       DO IELEM = 1,NELEM
00316 !
00317         K1 = -PHIEL(IELEM,1)
00318         K2 = -PHIEL(IELEM,2)
00319         K3 = -PHIEL(IELEM,3)
00320 !
00321 !       STARTS WITH N-SCHEME (EQUIVALENT TO LEO POSTMA'S IMPLEMENTATION)
00322 !       FIJ HERE LIKE LAMBDA(I,J) IN BOOK
00323 !
00324         F12=MAX(MIN(K1,-K2),0.D0)
00325         F23=MAX(MIN(K2,-K3),0.D0)
00326         F31=MAX(MIN(K3,-K1),0.D0)
00327         F21=MAX(MIN(K2,-K1),0.D0)
00328         F32=MAX(MIN(K3,-K2),0.D0)
00329         F13=MAX(MIN(K1,-K3),0.D0)
00330 !
00331         I1=IKLE(IELEM,1)
00332         I2=IKLE(IELEM,2)
00333         I3=IKLE(IELEM,3)
00334         FN1=FN%R(I1)
00335         FN2=FN%R(I2)
00336         FN3=FN%R(I3)
00337 !
00338         BETA1FI=F12*(FN1-FN2)+F13*(FN1-FN3)
00339         BETA2FI=F21*(FN2-FN1)+F23*(FN2-FN3)
00340         BETA3FI=F31*(FN3-FN1)+F32*(FN3-FN2)
00341 !
00342         FI=BETA1FI+BETA2FI+BETA3FI
00343 !
00344 !       NOW PSI-SCHEME
00345 !
00346 !       WHAT FOLLOWS IS INSPIRED FROM SUBROUTINE VC08AA
00347 !       WHERE FIJ IS LIJ
00348 !
00349 !       THIS LINE CHANGES THE SCHEME INTO N-SCHEME
00350 !       GO TO 1000
00351 !
00352         IF(FI.GT.0.D0) THEN
00353           IF(BETA1FI.GT.FI) THEN
00354             F12=F12*FI/BETA1FI
00355             F13=F13*FI/BETA1FI
00356           ELSEIF(BETA1FI.LT.0.D0) THEN
00357             F12=0.D0
00358             F13=0.D0
00359           ENDIF
00360           IF(BETA2FI.GT.FI) THEN
00361             F21=F21*FI/BETA2FI
00362             F23=F23*FI/BETA2FI
00363           ELSEIF(BETA2FI.LT.0.D0) THEN
00364             F21=0.D0
00365             F23=0.D0
00366           ENDIF
00367           IF(BETA3FI.GT.FI) THEN
00368             F31=F31*FI/BETA3FI
00369             F32=F32*FI/BETA3FI
00370           ELSEIF(BETA3FI.LT.0.D0) THEN
00371             F31=0.D0
00372             F32=0.D0
00373           ENDIF
00374         ELSEIF(FI.LT.0.D0) THEN
00375           IF(BETA1FI.LT.FI) THEN
00376             F12=F12*FI/BETA1FI
00377             F13=F13*FI/BETA1FI
00378           ELSEIF(BETA1FI.GT.0.D0) THEN
00379             F12=0.D0
00380             F13=0.D0
00381           ENDIF
00382           IF(BETA2FI.LT.FI) THEN
00383             F21=F21*FI/BETA2FI
00384             F23=F23*FI/BETA2FI
00385           ELSEIF(BETA2FI.GT.0.D0) THEN
00386             F21=0.D0
00387             F23=0.D0
00388           ENDIF
00389           IF(BETA3FI.LT.FI) THEN
00390             F31=F31*FI/BETA3FI
00391             F32=F32*FI/BETA3FI
00392           ELSEIF(BETA3FI.GT.0.D0) THEN
00393             F31=0.D0
00394             F32=0.D0
00395           ENDIF
00396         ELSE
00397           F12=0.D0
00398           F23=0.D0
00399           F31=0.D0
00400           F21=0.D0
00401           F32=0.D0
00402           F13=0.D0
00403         ENDIF
00404 !
00405 !1000    CONTINUE
00406 !
00407 !       ASSEMBLES FLUXES
00408 !       A DIFFERENCE WITH STANDARD DISTRIBUTIVE SCHEMES
00409 !
00410 !       HERE BEWARE CONVENTION ON FLOW
00411 !
00412 !       SEGMENT 1
00413         ISEG  = ELTSEG(IELEM,1)
00414         IF(ORISEG(IELEM,1).EQ.1) THEN
00415           FLOW(ISEG) = FLOW(ISEG) + F21 - F12
00416         ELSE
00417           FLOW(ISEG) = FLOW(ISEG) - F21 + F12
00418         ENDIF
00419 !       SEGMENT 2
00420         ISEG  = ELTSEG(IELEM,2)
00421         IF(ORISEG(IELEM,2).EQ.1) THEN
00422           FLOW(ISEG) = FLOW(ISEG) + F32 - F23
00423         ELSE
00424           FLOW(ISEG) = FLOW(ISEG) - F32 + F23
00425         ENDIF
00426 !       SEGMENT 3
00427         ISEG  = ELTSEG(IELEM,3)
00428         IF(ORISEG(IELEM,3).EQ.1) THEN
00429           FLOW(ISEG) = FLOW(ISEG) + F13 - F31
00430         ELSE
00431           FLOW(ISEG) = FLOW(ISEG) - F13 + F31
00432         ENDIF
00433       ENDDO
00434 !
00435 !-----------------------------------------------------------------------
00436 !
00437       ELSE
00438 !
00439         IF(LNG.EQ.1) THEN
00440           WRITE(LU,*) 'FLUX_EF_VF :'
00441           WRITE(LU,*) 'OPTION INCONNUE : ',IOPT
00442 !         IF(IOPT.EQ.3.AND..NOT.PRESENT(FN)) THEN
00443           IF(IOPT.EQ.3) THEN
00444             WRITE(LU,*) 'OPTION 3 : FONCTION CONVECTEE REQUISE'
00445           ENDIF
00446         ENDIF
00447         IF(LNG.EQ.2) THEN
00448           WRITE(LU,*) 'FLUX_EF_VF:'
00449           WRITE(LU,*) 'UNKNOWN OPTION: ',IOPT
00450 !         IF(IOPT.EQ.3.AND..NOT.PRESENT(FN)) THEN
00451           IF(IOPT.EQ.3) THEN
00452             WRITE(LU,*) 'OPTION 3: ADVECTED FUNCTION REQUIRED'
00453           ENDIF
00454         ENDIF
00455         CALL PLANTE(1)
00456         STOP
00457 !
00458       ENDIF
00459 !
00460 !-----------------------------------------------------------------------
00461 !
00462       RETURN
00463       END

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