prediv.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\telemac3d\prediv.f
00002 !
00146                      SUBROUTINE PREDIV
00147 !                    *****************
00148 !
00149      & ( PD, UP, VP, WP, INFO , BC , OPT, DIRSUR, DIRBOT, DIRLAT)
00150 !
00151 !***********************************************************************
00152 ! TELEMAC3D   V6P2                                   21/08/2010
00153 !***********************************************************************
00154 !
00155 !
00156 !
00157 !
00158 !
00159 !
00160 !
00161 !
00162 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00163 !| BC             |-->| LOGICAL, IF YES, BOUNDARY CONDITIONS
00164 !|                |   | ARE TAKEN INTO ACCOUNT (SEE CALL FLUPRI...)
00165 !|                |   | AND DIRICHLET VALUES ARE CONSIDERED
00166 !|                |   | IF NO, DIRICHLET VALUES ARE SET TO 0
00167 !|                |   | AND BOUNDARY CONDITIONS ARE DISCARDED
00168 !| DIRBOT         |-->| LOGICAL FOR DIRICHLET VALUES ON THE BOTTOM
00169 !| DIRLAT         |-->| LOGICAL FOR DIRICHLET VALUES ON LATERAL BOUNDARIES
00170 !| DIRSUR         |-->| LOGICAL FOR DIRICHLET VALUES AT THE SURFACE
00171 !| INFO           |-->| INFORMATION ASKED ON SOLVER ITERATIONS
00172 !| OPT            |-->| OPTION 1: DIVERGENCE IN THE REAL MESH
00173 !|                |   | 2: DIVERGENCE OBTAINED BY SUM OF FLUXES
00174 !|                |   | 3: DIVERGENCE OBTAINED BY SUM OF FLUXES
00175 !|                |   | AND LAPLACIAN IN THE TRANSFORMED MESH
00176 !|                |   | THAT WILL ALLOW A SPLITTING BETWEEN
00177 !|                |   | VERTICAL FLUXES AND FLUXES TANGENT
00178 !|                |   | TO PLANES
00179 !| PD             |<->| DYNAMIC PRESSURE
00180 !| UP             |-->| INTERMEDIATE VELOCITY FIELD
00181 !| VP             |-->| INTERMEDIATE VELOCITY FIELD
00182 !| WP             |-->| INTERMEDIATE VELOCITY FIELD
00183 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00184 !
00185       USE BIEF
00186       USE DECLARATIONS_TELEMAC
00187       USE DECLARATIONS_TELEMAC3D
00188 !
00189       IMPLICIT NONE
00190       INTEGER LNG,LU
00191       COMMON/INFO/LNG,LU
00192 !
00193 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00194 !
00195       INTEGER, INTENT(IN)             :: OPT
00196       TYPE(BIEF_OBJ), INTENT(INOUT)   :: PD
00197       TYPE(BIEF_OBJ), INTENT(IN)      :: UP, VP, WP
00198       LOGICAL, INTENT(IN)             :: INFO,BC,DIRSUR,DIRBOT,DIRLAT
00199 !
00200 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00201 !
00202       INTEGER IPOIN3,IPOIN2,IPTFR3,IPLAN,IS,IIS,SIZD,SIZX,ELMD
00203       CHARACTER(LEN=16) FORMUL
00204 !
00205 !=======================================================================
00206 ! RIGHT HAND SIDE VECTOR SEM3D = - DIV (INTERMEDIATE VELOCITY)
00207 !=======================================================================
00208 !
00209       IF(OPT.EQ.1) THEN
00210 !                     2 FILTERING CRUSHED ELEMENTS
00211         FORMUL='GRADF 2         '
00212         CALL VECTOR
00213      &   (SEM3D, '=',FORMUL(1:15)//'X', IELM3,-1.D0,UP, SVIDE,SVIDE,
00214      &    SVIDE, SVIDE, SVIDE, MESH3D, MSK, MASKEL)
00215         CALL VECTOR
00216      &   (SEM3D, '+',FORMUL(1:15)//'Y', IELM3,-1.D0,VP, SVIDE,SVIDE,
00217      &    SVIDE, SVIDE, SVIDE, MESH3D, MSK, MASKEL)
00218         CALL VECTOR
00219      &   (SEM3D, '+',FORMUL(1:15)//'Z', IELM3,-1.D0,WP, SVIDE,SVIDE,
00220      &    SVIDE, SVIDE, SVIDE, MESH3D, MSK, MASKEL)
00221 !
00222         IF(BC) THEN
00223           CALL FLUPRI(SEM3D%R,1.D0,UP%R,VP%R,WP%R,
00224      &                X,Y,Z,IKLE3%I,MESH3D%NELEM,
00225      &                MESH3D%NELMAX,NELEM2,NPOIN2,NPOIN3,
00226      &                MESH2D%W%R(         1:  NELEM2),
00227      &                MESH2D%W%R(  NELEM2+1:2*NELEM2),
00228      &                MESH2D%W%R(2*NELEM2+1:3*NELEM2)  )
00229 !         TERME EN D(ZS)/DT (MAIS ATTENTION DSSUDT=(H-HN)/DT)
00230           CALL VECTOR(T2_01,'=', 'MASVEC          ',IELM2H,-1.D0,DSSUDT,
00231      &                SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH2D,MSK,MASKEL)
00232           CALL OV('X=X+Y   ',SEM3D%R(NPOIN3-NPOIN2+1:NPOIN3),
00233      &                       T2_01%R,T2_01%R,0.D0,NPOIN2)
00234 !         MANQUE LE TERME EN D(ZF)/DT
00235         ENDIF
00236 !
00237 !       SOURCES INSIDE THE DOMAIN : DIV(U) WILL NOT BE 0
00238 !       WITH OPTIONS 2 AND 3, THIS IS ALREADY IN RHS_PRESSURE
00239 !
00240         IF(NSCE.GE.1) THEN
00241           DO IS=1,NSCE
00242             IIS=IS
00243 !           HERE SOURCES IN THE NON ASSEMBLED (PARCOM) FORM
00244 !           SEE SOURCES_SINKS
00245             IF(NCSIZE.GT.1) IIS=IS+NSCE
00246             CALL OS('X=X+Y   ',X=SEM3D,Y=SOURCES%ADR(IIS)%P)
00247           ENDDO
00248         ENDIF
00249 !
00250       ELSEIF(OPT.EQ.2.OR.OPT.EQ.3) THEN
00251 !
00252         CALL RHS_PRESSURE(SEM3D,UP,VP,WP,IELM3,DM1,GRAZCO,SVIDE,
00253      &                    MESH3D,MSK,MASKEL,FLUEXT,NSCE,RAIN,PLUIE,
00254      &                    SOURCES,GRADZF,VOLU2D,DSSUDT,
00255      &                    NPOIN2,NPOIN3,NPLAN)
00256 !
00257       ELSE
00258 !
00259         IF(LNG.EQ.1) WRITE(LU,*) 'OPTION INCONNUE DANS PREDIV: ',OPT
00260         IF(LNG.EQ.2) WRITE(LU,*) 'UNKNOWN OPTION IN PREDIV: ',OPT
00261         CALL PLANTE(1)
00262         STOP
00263 !
00264       ENDIF
00265 !
00266 !=======================================================================
00267 ! DIFFUSION MATRIX AND BOUNDARY TERMS (BC)
00268 !=======================================================================
00269 !
00270       CALL CPSTVC(SEM3D,T3_01)
00271       CALL OS('X=C     ', X=T3_01,C=1.D0)
00272 !
00273       IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00274         FORMUL='MATDIF       MON'
00275       ELSE
00276         FORMUL='MATDIF          '
00277       ENDIF
00278 !
00279       IF(OPT.EQ.3) THEN
00280         FORMUL='MATDIF*  1234   '
00281         IF(SIGMAG.OR.OPTBAN.EQ.1) FORMUL(7:7)='2'
00282       ENDIF
00283 !
00284       CALL MATRIX(MDIFF,'M=N     ',FORMUL,IELM3,IELM3,1.D0,
00285      &            T3_01,T3_01,T3_01,
00286      &            SVIDE,SVIDE,SVIDE,MESH3D,MSK,MASKEL)
00287 !
00288 !=======================================================================
00289 ! IMPOSED (DIRICHLET) BC + PRECONDITIONING + LINEAR EQUATION SOLVER
00290 !=======================================================================
00291 !
00292 ! IMPOSED BC
00293 !
00294       DO IPOIN3 = 1,NPOIN3
00295         IT1%I(IPOIN3) = KDDL
00296         IT2%I(IPOIN3) = 0
00297         T3_03%R(IPOIN3) = 0.D0
00298       ENDDO
00299 !
00300 !     BORDS LATERAUX
00301 !
00302       IF(DIRLAT) THEN
00303       IF(BC) THEN
00304         DO IPTFR3 = 1,NPTFR3
00305           IF(LIPBOL%I(IPTFR3).EQ.KENT.OR.LIPBOL%I(IPTFR3).EQ.KADH) THEN
00306             IT1%I(NBOR3%I(IPTFR3)) = KDIR
00307             T3_03%R(NBOR3%I(IPTFR3)) = PBORL%R(IPTFR3)
00308           ENDIF
00309         ENDDO
00310       ELSE
00311 !       PBORL REPLACED BY 0
00312         DO IPTFR3 = 1,NPTFR3
00313           IF(LIPBOL%I(IPTFR3).EQ.KENT.OR.LIPBOL%I(IPTFR3).EQ.KADH) THEN
00314             IT1%I(NBOR3%I(IPTFR3)) = KDIR
00315             T3_03%R(NBOR3%I(IPTFR3)) = 0.D0
00316           ENDIF
00317         ENDDO
00318       ENDIF
00319       ENDIF
00320 !
00321 !     BOTTOM
00322 !
00323       IF(DIRBOT) THEN
00324       IF(BC) THEN
00325         DO IPOIN2 = 1,NPOIN2
00326           IF(LIPBOF%I(IPOIN2).EQ.KENT.OR.LIPBOF%I(IPOIN2).EQ.KADH) THEN
00327             IT1%I(IPOIN2) = KDIR
00328             T3_03%R(IPOIN2) = PBORF%R(IPOIN2)
00329           ENDIF
00330         ENDDO
00331       ELSE
00332         DO IPOIN2 = 1,NPOIN2
00333           IF(LIPBOF%I(IPOIN2).EQ.KENT.OR.LIPBOF%I(IPOIN2).EQ.KADH) THEN
00334             IT1%I(IPOIN2) = KDIR
00335             T3_03%R(IPOIN2) = 0.D0
00336           ENDIF
00337         ENDDO
00338       ENDIF
00339       ENDIF
00340 !
00341 !     FREE SURFACE
00342 !
00343       IF(DIRSUR) THEN
00344       IF(NCSIZE.GT.1) THEN
00345         DO IPOIN2 = 1,NPOIN2
00346           IF(LIPBOS%I(IPOIN2).EQ.KENT.OR.LIPBOS%I(IPOIN2).EQ.KADH) THEN
00347             IPOIN3=NPOIN3-NPOIN2+IPOIN2
00348             IT1%I(IPOIN3) = KDIR
00349             T3_03%R(IPOIN3) = PBORS%R(IPOIN2)
00350             MDIFF%D%R(IPOIN3)=MESH3D%FAC%R(IPOIN3)
00351           ENDIF
00352         ENDDO
00353       ELSE
00354         DO IPOIN2 = 1,NPOIN2
00355           IF(LIPBOS%I(IPOIN2).EQ.KENT.OR.LIPBOS%I(IPOIN2).EQ.KADH) THEN
00356             IPOIN3=NPOIN3-NPOIN2+IPOIN2
00357             IT1%I(IPOIN3) = KDIR
00358             T3_03%R(IPOIN3) = PBORS%R(IPOIN2)
00359             MDIFF%D%R(IPOIN3)=1.D0
00360           ENDIF
00361         ENDDO
00362       ENDIF
00363       ENDIF
00364 !
00365 !-----------------------------------------------------------------------
00366 !
00367 !     CRUSHED ELEMENTS AND TIDAL FLATS: FIRST TREATED AS DIRICHLET
00368 !     A CRUSHED POINT HAS NO DIAGONAL
00369 !
00370       IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00371         IF(NCSIZE.GT.1) THEN
00372           CALL OS('X=Y     ',X=T3_02,Y=MDIFF%D)
00373 !         T3_02 WILL BE THE ASSEMBLED DIAGONAL
00374           CALL PARCOM(T3_02,2,MESH3D)
00375           DO IPOIN3=1,NPOIN3
00376             IF(T3_02%R(IPOIN3).LT.1.D-10) THEN
00377               MDIFF%D%R(IPOIN3)= MESH3D%FAC%R(IPOIN3)
00378               IT1%I(IPOIN3)    = KDIR
00379               IT2%I(IPOIN3)    = 1
00380               T3_03%R(IPOIN3)  = PD%R(IPOIN3)
00381             ENDIF
00382           ENDDO
00383         ELSE
00384           DO IPOIN3=1,NPOIN3
00385             IF(MDIFF%D%R(IPOIN3).LT.1.D-10) THEN
00386               MDIFF%D%R(IPOIN3)= 1.D0
00387               IT1%I(IPOIN3)    = KDIR
00388               IT2%I(IPOIN3)    = 1
00389               T3_03%R(IPOIN3)  = PD%R(IPOIN3)
00390             ENDIF
00391           ENDDO
00392         ENDIF
00393       ENDIF
00394 !
00395 !-----------------------------------------------------------------------
00396 !
00397       CALL DIRI01(PD,MDIFF,SEM3D,T3_03,IT1%I,T3_01,T3_02,MESH3D,
00398      &            KDIR,MSK,MASKPT)
00399 !
00400 !     SOLVE THE LINEAR EQUATION SYSTEM MDIFF * PD = SEM3D
00401 !
00402       IF(NPLAN.EQ.2.AND.MDIFF%STO.EQ.3.AND.
00403      &   MDIFF%ELMLIN.EQ.41.AND.MDIFF%ELMCOL.EQ.41) THEN
00404 !
00405 !       AS THE UPPER PLANE IS DIRICHLET, THIS IS A 2D PROBLEM
00406 !
00407 !       SEM3D AND PD TURNED INTO 2D
00408         CALL CPSTVC(U2D,PD)
00409         CALL CPSTVC(U2D,SEM3D)
00410 !
00411 !       PROPERTIES OF MATRIX MODIFIED
00412 !       THE FIRST SEGMENTS IN A 3D MATRIX ARE THE HORIZONTAL
00413 !       SEGMENTS OF THE BOTTOM, HENCE JUST LIKE A 2D MATRIX
00414 !       THE FIRST POINTS IN THE DIAGONAL ARE THE BOTTOM POINTS
00415 !
00416         SIZD=MDIFF%D%DIM1
00417         SIZX=MDIFF%X%DIM1
00418         ELMD=MDIFF%D%ELM
00419         MDIFF%D%DIM1=NPOIN2
00420         MDIFF%X%DIM1=MESH2D%NELEM
00421         MDIFF%ELMLIN=11
00422         MDIFF%ELMCOL=11
00423         MDIFF%D%ELM=11
00424 !       PRECONDITIONING 17 WILL NOT BE ACCEPTED IN 2D
00425         IS=SLVPOI%PRECON
00426         IF(17*(IS/17).EQ.IS) SLVPOI%PRECON=IS/17
00427 !
00428         CALL SOLVE(PD,MDIFF,SEM3D,TRAV2,SLVPOI,INFO,MESH2D,
00429      &             MAT2D%ADR(1)%P)
00430 !
00431 !       RESTORING POSSIBLE PRECONDITIONING 17
00432         SLVPOI%PRECON=IS
00433 !
00434 !       OLD PROPERTIES RESTORED
00435         MDIFF%ELMLIN=41
00436         MDIFF%ELMCOL=41
00437         MDIFF%D%DIM1=SIZD
00438         MDIFF%X%DIM1=SIZX
00439         MDIFF%D%ELM=ELMD
00440         CALL CPSTVC(U,PD)
00441         CALL CPSTVC(U,SEM3D)
00442 !
00443       ELSE
00444 !
00445 !       REAL 3D PROBLEM
00446 !
00447         CALL SOLVE(PD,MDIFF,SEM3D,TRAV3,SLVPOI,INFO,MESH3D,MTRA2)
00448 !
00449       ENDIF
00450 !
00451 !=======================================================================
00452 !
00453 !   POINTS WITH NO DIAGONAL ARE EQUAL TO THE POINT ABOVE
00454 !   THEY HAD BEEN TREATED ABOVE AS DIRICHLET.
00455 !
00456 !=======================================================================
00457 !
00458       IF(SIGMAG.OR.OPTBAN.EQ.1) THEN
00459         DO IPLAN=NPLAN-1,1,-1
00460           DO IPOIN2=1,NPOIN2
00461             IPOIN3=IPOIN2+(IPLAN-1)*NPOIN2
00462             IF(IT2%I(IPOIN3).EQ.1) THEN
00463 !             VALUE OF THE UPPER POINT IS COPIED
00464               PD%R(IPOIN3)=PD%R(IPOIN3+NPOIN2)
00465             ENDIF
00466           ENDDO
00467         ENDDO
00468       ENDIF
00469 !
00470 !-----------------------------------------------------------------------
00471 !
00472       RETURN
00473       END

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