pre322
C PRE322 SOURCE OF166741 24/10/03 21:15:34 12022 & ICEN,IFACE,IFACEL,INORM, & IROC, IGRROC, IALROC, & IVITC, IGRVC, IALVC, & IPC ,IGRPC, IALPC, & IYC ,IGRYC, IALYC, & ISCAC ,IGRSC, IALSC, & IGAMC, & DELTAT, & IROF,IVITF,IPF,IYF,ISCAF, & LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : PRE322 C C DESCRIPTION : Voir PRE32 C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF C C************************************************************************ C C C APPELES (Outils) : KRIPAD, LICHT C C APPELES (Calcul) : AUCUN C C C************************************************************************ C C ENTREES C C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps C sinon 1er ordre en temps C LGAMC,LYC,LSCAC : LOGICAL : si .TRUE. IGAMC, IYC, ISCAC sont C pointeurs de CHPOINTS C C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE C C ICEN : MELEME de 'POI1' SPG des CENTRES C C IFACE : MELEME de 'POI1' SPG des FACES C C IFACEL : MELEME de 'SEG3' avec C CENTRE d'Elt "gauche" C CENTRE de Face C CENTRE d'Elt "droite" C C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2) C C INORM : CHPOINT des cosinus directeurs de normales aux faces C C 2) Autres pointeurs C C IROC : CHPOINT "CENTRE" contenant la masse volumique RHO C C IGRROC : CHPOINT "CENTRE" contenant le gradient de la C masse volumique RHO (2 composantes) C C IALROC : CHPOINT "CENTRE" contenant le limiteur du gradient C de la masse volumique C C IVITC : CHPOINT "CENTRE" contenant la vitesse UX, UY ; C C IGRVC : CHPOINT "CENTRE" contenant le gradient de la C vitesse (4 composantes) C C IALVC : CHPOINT "CENTRE" contenant le limiteur du gradient C de la vitesse (2 composantes) C C IPC : CHPOINT "CENTRE" contenat la pression P; C C IGRPC : CHPOINT "CENTRE" contenant le gradient de la C pression (2 composantes) C C IALPC : CHPOINT "CENTRE" contenant le limiteur du gradient C de la pression C C IYC : CHPOINT "CENTRE" contenat les fractions massiques ; C C IGRYC : CHPOINT "CENTRE" contenant les gradient des fr.mass. ; C C IALPC : CHPOINT "CENTRE" contenant les limiteurs des gradients C des fr.mass. ; C C ISCAC : CHPOINT "CENTRE" contenat les scalaires passifs ; C C IGRSC : CHPOINT "CENTRE" contenant les gradient des scalaires C passifs ; C C IALPC : CHPOINT "CENTRE" contenant les limiteurs des gradients C des scalaires passifs ; C C IGAMC : CHPOINT "CENTRE" contenat le "Gamma" du gaz ; C C 3) C C DELTAT : REAL*8, encrement en temps pour calculer la prediction C C C SORTIES C C C IROF : MCHAML defini sur le MELEME de pointeur IFACEL, C contenant la masse volumique RHO C C IVITF : MCHAML defini sur le MELEME de pointeur IFACEL, C contenant la vitesse UN, UT dans le repaire local C (n,t) et defini sur le MELEME de pointeur IFACE, C contenant les cosinus directeurs du repere local C C IPF : MCHAML defini sur le MELEME de pointeur IFACEL, C contenant la pression P C C IYF : MCHAML defini sur le MELEME de pointeur IFACEL, C contenant les fractions massiques; C C ISCAF : MCHAML defini sur le MELEME de pointeur IFACEL, C contenant les scalaires passifs; C C LOGAN : anomalie detectee C C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité C negative a été detectée -> en interactif le C programme s'arrete en GIBIANE C (erreur stocké en MESERR et VALER) C C LOGBOR : (LOGICAL): si .TRUE. un y a ete detecte C dehor 1 et 3 (sa valeur stockée en MESERR et VALER; C en VAL1 et en VAL2 on stocke 1.0 et 3.0) C C MESERR C VALER C VAL1, C VAL2 : pour les messages d'erreur C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : Créée le 21.12.98. C C 17.02.2000: transport des scalaires passifs C C************************************************************************ C C C ATTENTION: Cet programme marche si le MAILLAGE est convex; C si non il faut changer l'algoritme de calcul de C l'orientation des normales aux faces. C C La positivité n'est pas controlle parce que c'est déjà fait C dans l'operateur PRIM C C C************************************************************************ C IMPLICIT INTEGER(I-N) C**** Les variables C INTEGER ICEN, IFACE, IFACEL, INORM,IROC, IGRROC, IALROC & , IVITC, IGRVC, IALVC & , IPC ,IGRPC, IALPC & , IYC, IGRYC, IALYC & , ISCAC, IGRSC, IALSC & , IGAMC & , IROF, IVITF, IPF, IYF, ISCAF & , IGEOM, NFAC, NCEN & , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE & , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1 & , IDIMP1, INDCEL, I1, NESP, NSCA, NMA REAL*8 VALER, VAL1, VAL2, XG, YG, ZG, XC, YC, ZC, XD, YD, ZD, & DELTAT & ,DXG, DYG, DZG, DXD, DYD, DZD,ORIENT & , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ & , ROG, PG, GAMG, UXG, UYG, UZG, UNG, UTG, UVG & , ROD, PD, UXD, UYD, UZD, UND, UTD, UVD & , VALCEL, DCEL, ALCEL & , DROX, DROY, DROZ, DUXX, DUXY, DUXZ, DUYX, DUYY, & DUYZ, DUZX, DUZY, DUZZ, DPX, DPY, DPZ & , DRO, DUX, DUY, DUZ, DP & , ALPHA, DYMAS, SUMY CHARACTER*(40) MESERR CHARACTER*(8) TYPE LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGI2 & , LGAMC,LYC,LSCAC C C**** Les Includes C -INC SMCOORD -INC PPARAM -INC CCOPTIO -INC SMCHPOI POINTEUR MPROC.MPOVAL, MPGRR.MPOVAL, & MPVITC.MPOVAL, MPGRV.MPOVAL, & MPPC.MPOVAL, MPGRP.MPOVAL, & MPYC.MPOVAL, MPGRY.MPOVAL, & MPSCAC.MPOVAL, MPGRS.MPOVAL, & MPGAMC.MPOVAL, MPNORM.MPOVAL, & MPROP.MPOVAL, MPPP.MPOVAL, MPVITP.MPOVAL, & MPYP.MPOVAL, MPSCAP.MPOVAL -INC SMCHAML POINTEUR MCHAMY.MCHAML, MCHAMS.MCHAML POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL, & MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL, & MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL POINTEUR MELVUN.MELVAL, MELVUT.MELVAL, MELVUV.MELVAL POINTEUR MELRO.MELVAL, MELP.MELVAL POINTEUR MELALR.MPOVAL, & MELALV.MPOVAL, & MELALP.MPOVAL, & MELALY.MPOVAL, & MELALS.MPOVAL -INC SMLENTI -INC SMELEME C C**** Segments des fractions massiques gauche et droit C SEGMENT FRAMAS REAL*8 FRAMG(NMA), FRAMD(NMA) ENDSEGMENT POINTEUR SCALPA.FRAMAS C C**** Initialisation des parametres d'erreur déjà faite, i.e. C C LOGNEG = .FALSE. C LOGBOR = .FALSE. C MESERR = ' ' C MOTERR(1:40) = MESERR(1:40) C VALER = 0.0D0 C VAL1 = 0.0D0 C VAL2 = 0.0D0 C C C**** KRIPAD pour la correspondance global/local de centre C C C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements C C Si i est le numero global d'un noeud de ICEN, C MLENT1.LECT(i) contient sa position, i.e. C C I = numero global du noeud centre C MLENT1.LECT(i) = numero local du noeud centre C C MLENT1 déjà activé, i.e. C C SEGACT MLENT1 C C**** Activation de CHPOINTs C C densité + grad + limiteur C vitesse + grad + limiteur C pression + grad + limiteur C fract.mass + grad + limiteur C gamma C cosinus directeurs des normales aux surface C C C**** Les MPOVALs 'Prediction' C IF(LOGTEM)THEN SEGINI, MPROP = MPROC SEGINI, MPPP = MPPC SEGINI, MPVITP = MPVITC ELSE MPROP = MPROC MPPP = MPPC MPVITP = MPVITC ENDIF C C**** Les Limiteurs C C C C**** Les MPOVAL sont déjà activés i.e.: C C SEGACT MPROC C SEGACT MPGRR C SEGACT MPIALR C SEGACT MPVITC C SEGACT MPGRV C SEGACT MPIALV C SEGACT MPPC C SEGACT MPGRP C SEGACT MPIALP C SEGACT MPGAMC C SEGACT MPNORM C C**** Le MELEME FACEL C IPT1 = IFACEL IPT2 = IFACE SEGACT IPT1 SEGACT IPT2 NFAC = IPT1.NUM(/2) C C**** Creation de MCHAMLs contenant les etats gauche et droite, C C i.e.: C C vitesse + cosinus directors du repere local C densité C pression C C**** Cosinus directors du repere local et vitesse C C Les cosinus directeurs C N1 = 2 N3 = 6 L1 = 28 SEGINI MCHEL1 IVITF = MCHEL1 MCHEL1.TITCHE = 'U ' MCHEL1.IMACHE(1) = IFACE MCHEL1.IMACHE(2) = IFACEL MCHEL1.CONCHE(1) = '(n,t,v)in(x,y,z)' MCHEL1.CONCHE(2) = ' U in (n,t,v) ' C C**** Valeurs des cosinus definies par respect au repair global, i.e. C MCHEL1.INFCHE(1,1) = 2 MCHEL1.INFCHE(1,3) = NIFOUR MCHEL1.INFCHE(1,4) = 0 MCHEL1.INFCHE(1,5) = 0 MCHEL1.INFCHE(1,6) = 1 MCHEL1.IFOCHE = IFOUR C C**** Valeurs de vitesse definies par respect au repair local, i.e. C MCHEL1.INFCHE(2,1) = 1 MCHEL1.INFCHE(2,3) = NIFOUR MCHEL1.INFCHE(2,4) = 0 MCHEL1.INFCHE(2,5) = 0 MCHEL1.INFCHE(2,6) = 1 C C**** Le cosinus directeurs C N1PTEL = 1 N1EL = NFAC N2PTEL = 0 N2EL = 0 C C**** MCHAML a N2 composantes: C C cosinus directeurs du repere local (n,t1) C C IDIM = 3 -> 9 composantes C N2 = 9 SEGINI MCHAM1 MCHEL1.ICHAML(1) = MCHAM1 MCHAM1.NOMCHE(1) = 'NX ' MCHAM1.NOMCHE(2) = 'NY ' MCHAM1.NOMCHE(3) = 'NZ ' MCHAM1.NOMCHE(4) = 'TX ' MCHAM1.NOMCHE(5) = 'TY ' MCHAM1.NOMCHE(6) = 'TZ ' MCHAM1.NOMCHE(7) = 'VX ' MCHAM1.NOMCHE(8) = 'VY ' MCHAM1.NOMCHE(9) = 'VZ ' MCHAM1.TYPCHE(1) = 'REAL*8 ' MCHAM1.TYPCHE(2) = 'REAL*8 ' MCHAM1.TYPCHE(3) = 'REAL*8 ' MCHAM1.TYPCHE(4) = 'REAL*8 ' MCHAM1.TYPCHE(5) = 'REAL*8 ' MCHAM1.TYPCHE(6) = 'REAL*8 ' MCHAM1.TYPCHE(7) = 'REAL*8 ' MCHAM1.TYPCHE(8) = 'REAL*8 ' MCHAM1.TYPCHE(9) = 'REAL*8 ' SEGINI MELVNX SEGINI MELVNY SEGINI MELVNZ SEGINI MELT1X SEGINI MELT1Y SEGINI MELT1Z SEGINI MELT2X SEGINI MELT2Y SEGINI MELT2Z MCHAM1.IELVAL(1) = MELVNX MCHAM1.IELVAL(2) = MELVNY MCHAM1.IELVAL(3) = MELVNZ MCHAM1.IELVAL(4) = MELT1X MCHAM1.IELVAL(5) = MELT1Y MCHAM1.IELVAL(6) = MELT1Z MCHAM1.IELVAL(7) = MELT2X MCHAM1.IELVAL(8) = MELT2Y MCHAM1.IELVAL(9) = MELT2Z SEGDES MCHAM1 C C**** Vitesse C N1EL = NFAC N1PTEL = 3 N2EL = 0 N2PTEL = 0 C C**** MCHAML a N2 composantes: C C IDIM = 3 -> 3 composantes C N2 = 3 SEGINI MCHAM1 MCHEL1.ICHAML(2) = MCHAM1 SEGDES MCHEL1 MCHAM1.NOMCHE(1) = 'UN ' MCHAM1.NOMCHE(2) = 'UT ' MCHAM1.NOMCHE(3) = 'UV ' MCHAM1.TYPCHE(1) = 'REAL*8 ' MCHAM1.TYPCHE(2) = 'REAL*8 ' MCHAM1.TYPCHE(3) = 'REAL*8 ' SEGINI MELVUN SEGINI MELVUT SEGINI MELVUV MCHAM1.IELVAL(1) = MELVUN MCHAM1.IELVAL(2) = MELVUT MCHAM1.IELVAL(3) = MELVUV SEGDES MCHAM1 C C**** Densite C N1 = 1 N3 = 6 L1 = 15 SEGINI MCHEL2 IROF = MCHEL2 MCHEL2.IMACHE(1) = IFACEL MCHEL2.TITCHE = 'RO ' MCHEL2.CONCHE(1) = ' ' C C**** Valeurs independente du repére, i.e. C MCHEL2.INFCHE(1,1) = 0 MCHEL2.INFCHE(1,3) = NIFOUR MCHEL2.INFCHE(1,4) = 0 MCHEL2.INFCHE(1,5) = 0 MCHEL2.INFCHE(1,6) = 1 MCHEL2.IFOCHE = IFOUR N2 = 1 SEGINI MCHAM1 MCHEL2.ICHAML(1) = MCHAM1 SEGDES MCHEL2 MCHAM1.NOMCHE(1) = 'SCAL ' MCHAM1.TYPCHE(1) = 'REAL*8 ' SEGINI MELRO MCHAM1.IELVAL(1) = MELRO SEGDES MCHAM1 C C**** Pression C MCHEL1 = IROF SEGINI, MCHEL2 = MCHEL1 IPF = MCHEL2 MCHEL2.TITCHE = 'P ' C C**** MCHAM1 = MCHAML de la densite C SEGINI, MCHAM2 = MCHAM1 MCHEL2.ICHAML(1) = MCHAM2 SEGDES MCHEL2 SEGINI MELP MCHAM2.IELVAL(1) = MELP SEGDES MCHAM2 C C**** Le CHAMELEM des fractions massiques C IF(LYC)THEN MCHPO1 = IYC SEGACT MCHPO1 MSOUP1 = MCHPO1.IPCHP(1) SEGDES MCHPO1 SEGACT MSOUP1 NESP = MSOUP1.NOCOMP(/2) MPYC = MSOUP1.IPOVAL SEGACT MPYC NESP = MPYC.VPOCHA(/2) MCHEL1 = IROF SEGINI, MCHEL2 = MCHEL1 IYF = MCHEL2 MCHEL2.TITCHE = 'Y ' N2 = NESP SEGINI MCHAMY MCHEL2.ICHAML(1) = MCHAMY SEGDES MCHEL2 N1EL = NFAC N1PTEL = 3 N2EL = 0 N2PTEL = 0 DO I1 = 1, NESP SEGINI MELVA1 MCHAMY.IELVAL(I1) = MELVA1 MCHAMY.NOMCHE(I1) = MSOUP1.NOCOMP(I1) MCHAMY.TYPCHE(I1) = 'REAL*8 ' ENDDO SEGDES MSOUP1 C NMA = NESP SEGINI FRAMAS IF(LOGTEM)THEN SEGINI, MPYP = MPYC ELSE MPYP = MPYC ENDIF ELSE NESP = 0 IYF = 0 ENDIF C C**** Le CHAMELEM des scalaires passifs C IF(LSCAC)THEN MCHPO1 = ISCAC SEGACT MCHPO1 MSOUP1 = MCHPO1.IPCHP(1) SEGDES MCHPO1 SEGACT MSOUP1 MPSCAC = MSOUP1.IPOVAL SEGACT MPSCAC NSCA = MPSCAC.VPOCHA(/2) MCHEL1 = IROF SEGINI, MCHEL2 = MCHEL1 ISCAF = MCHEL2 MCHEL2.TITCHE = 'SCALPASS ' N2 = NSCA SEGINI MCHAMS MCHEL2.ICHAML(1) = MCHAMS SEGDES MCHEL2 N1EL = NFAC N1PTEL = 3 N2EL = 0 N2PTEL = 0 DO I1 = 1, NSCA, 1 SEGINI MELVA1 MCHAMS.IELVAL(I1) = MELVA1 MCHAMS.NOMCHE(I1) = MSOUP1.NOCOMP(I1) MCHAMS.TYPCHE(I1) = 'REAL*8 ' ENDDO SEGDES MSOUP1 C NMA = NSCA SEGINI SCALPA IF(LOGTEM)THEN SEGINI, MPSCAP = MPSCAC ELSE MPSCAP = MPSCAC ENDIF ELSE NSCA = 0 ISCAF = 0 ENDIF C C**** Donc on a aussi actives le chpoints de fractions massiques C C SEGACT MPYC C SEGACT MPGRY C SEGACT MPIALY C C C SEGACT MPSCAC C SEGACT MPGRS C SEGACT MPIALS C C C*********************************************************************** C********* PREDICTION ************************************************** C*********************************************************************** C C**** Prediction avec gradients limités C C IF(LOGTEM)THEN C IPT3 = ICEN SEGACT IPT3 NCEN = IPT3.NUM(/2) DO NLCE = 1, NCEN ROG = MPROP.VPOCHA(NLCE,1) UXG = MPVITP.VPOCHA(NLCE,1) UYG = MPVITP.VPOCHA(NLCE,2) UZG = MPVITP.VPOCHA(NLCE,3) PG = MPPP.VPOCHA(NLCE,1) DROX = MPGRR.VPOCHA(NLCE,1)*MELALR.VPOCHA(NLCE,1) DROY = MPGRR.VPOCHA(NLCE,2)*MELALR.VPOCHA(NLCE,1) DROZ = MPGRR.VPOCHA(NLCE,3)*MELALR.VPOCHA(NLCE,1) DUXX = MPGRV.VPOCHA(NLCE,1)*MELALV.VPOCHA(NLCE,1) DUXY = MPGRV.VPOCHA(NLCE,2)*MELALV.VPOCHA(NLCE,1) DUXZ = MPGRV.VPOCHA(NLCE,3)*MELALV.VPOCHA(NLCE,1) DUYX = MPGRV.VPOCHA(NLCE,4)*MELALV.VPOCHA(NLCE,2) DUYY = MPGRV.VPOCHA(NLCE,5)*MELALV.VPOCHA(NLCE,2) DUYZ = MPGRV.VPOCHA(NLCE,6)*MELALV.VPOCHA(NLCE,2) DUZX = MPGRV.VPOCHA(NLCE,7)*MELALV.VPOCHA(NLCE,3) DUZY = MPGRV.VPOCHA(NLCE,8)*MELALV.VPOCHA(NLCE,3) DUZZ = MPGRV.VPOCHA(NLCE,9)*MELALV.VPOCHA(NLCE,3) DPX = MPGRP.VPOCHA(NLCE,1)*MELALP.VPOCHA(NLCE,1) DPY = MPGRP.VPOCHA(NLCE,2)*MELALP.VPOCHA(NLCE,1) DPZ = MPGRP.VPOCHA(NLCE,3)*MELALP.VPOCHA(NLCE,1) GAMG = MPGAMC.VPOCHA(NLCE,1) DRO = UXG * DROX + ROG * ( DUXX + DUYY + DUZZ) & + UYG * DROY + UZG * DROZ DUX = UXG * DUXX + DPX / ROG + UYG * DUXY & + UZG * DUXZ DUY = UXG * DUYX + UYG * DUYY + DPY / ROG & + UZG * DUYZ DUZ = UXG * DUZX + UYG * DUZY + UZG * DUZZ & + DPZ / ROG & + UXG * DPX + UYG * DPY + UZG * DPZ C MPROP.VPOCHA(NLCE,1) = ROG - DELTAT * DRO MPVITP.VPOCHA(NLCE,1) = UXG - DELTAT * DUX MPVITP.VPOCHA(NLCE,2) = UYG - DELTAT * DUY MPVITP.VPOCHA(NLCE,3) = UZG - DELTAT * DUZ DO I1 = 1, NESP INDCEL = 3 * (I1-1) + 1 MPYP.VPOCHA(NLCE,I1) = MPYC.VPOCHA(NLCE,I1) - & DELTAT * DYMAS ENDDO C DO I1 = 1, NSCA, 1 INDCEL = 3 * (I1 - 1) + 1 MPSCAP.VPOCHA(NLCE,I1) = MPSCAC.VPOCHA(NLCE,I1) - & DELTAT * DYMAS ENDDO ENDDO C ENDIF C C C*********************************************************************** C********* CORRECTION ************************************************** C*********************************************************************** C C**** Boucle sur le faces C IDIMP1 = IDIM + 1 DO NLCF = 1, NFAC C C******* NLCF = numero local du centre de face C NGCF = numero global du centre de face C NGCEG = numero global du centre ELT "gauche" C NLCEG = numero local du centre ELT "gauche" C NGCED = numero global du centre ELT "droite" C NLCED = numero local du centre ELT "droite" C NGCEG = IPT1.NUM(1,NLCF) NGCF = IPT1.NUM(2,NLCF) NGCED = IPT1.NUM(3,NLCF) NLCEG = MLENT1.LECT(NGCEG) NLCED = MLENT1.LECT(NGCED) C C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF) C NGCF1 = IPT2.NUM(1,NLCF) IF( NGCF1 .NE. NGCF) THEN LOGAN = .TRUE. MESERR(1:40) = 'PRET, subroutine pre321.eso ' GOTO 9999 ENDIF C C******* Cosinus directeurs des NORMALES aux faces C C On impose que les normales sont direct "Gauche" -> "Centre" C INDCEL = (NGCEG-1)*IDIMP1 XG = XCOOR(INDCEL+1) YG = XCOOR(INDCEL+2) ZG = XCOOR(INDCEL+3) INDCEL = (NGCF-1)*IDIMP1 XC = XCOOR(INDCEL + 1) YC = XCOOR(INDCEL + 2) ZC = XCOOR(INDCEL+3) INDCEL = (NGCED-1)*IDIMP1 XD = XCOOR(INDCEL+1) YD = XCOOR(INDCEL+2) ZD = XCOOR(INDCEL+3) DXG = XC - XG DYG = YC - YG DZG = ZC - ZG DXD = XC - XD DYD = YC - YD DZD = ZC - ZD C C******* On calcule le sign du pruduit scalare C (Normales de Castem) * (vecteur "gauche" -> "centre") C CNX = MPNORM.VPOCHA(NLCF,7) CNY = MPNORM.VPOCHA(NLCF,8) CNZ = MPNORM.VPOCHA(NLCF,9) LOGAN = .TRUE. MESERR(1:30)= & 'PRET , subroutine pre322.eso. ' GOTO 9999 ENDIF C C********** Cosinus directeurs de tangente 1 C C C********** Cosinus directeurs de tangente 2 C C C C******* Les autres MELVALs C C C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0 C y_i \in (0,1) C Si non il faut le faire, en utilisant LOGBOR, C LOGNEG, VALER, VAL1, VAL2 C C C C******* NGCEG = NGCED -> Mur C IF( NGCEG .EQ. NGCED)THEN C C********** Sur le mur on fait de reconstruction sur l'etat gauche C C C********** Etat gauche C VALCEL = MPROP.VPOCHA(NLCEG, 1) ALCEL = MELALR.VPOCHA(NLCEG, 1) DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG + & MPGRR.VPOCHA(NLCEG, 2)*DYG + & MPGRR.VPOCHA(NLCEG, 3)*DZG ROG = VALCEL + ALCEL * DCEL C VALCEL = MPPP.VPOCHA(NLCEG, 1) ALCEL = MELALP.VPOCHA(NLCEG, 1) DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG + & MPGRP.VPOCHA(NLCEG, 2)*DYG + & MPGRP.VPOCHA(NLCEG, 3)*DZG PG = VALCEL + ALCEL * DCEL C LOGI2 = .FALSE. SUMY = 0.0D0 DO I1 = 1, NESP, 1 INDCEL = 3 * (I1-1) + 1 VALCEL = MPYP.VPOCHA(NLCEG,I1) ALCEL = MELALY.VPOCHA(NLCEG, I1) DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG + & MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG + & MPGRY.VPOCHA(NLCEG,INDCEL + 2 )*DZG ALCEL = VALCEL + ALCEL * DCEL SUMY = SUMY + ALCEL LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0) FRAMAS.FRAMG(I1) = ALCEL ENDDO LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0) C DO I1 = 1, NSCA, 1 INDCEL = 3 * (I1-1) + 1 VALCEL = MPSCAP.VPOCHA(NLCEG,I1) ALCEL = MELALS.VPOCHA(NLCEG, I1) DCEL = MPGRS.VPOCHA(NLCEG, INDCEL)*DXG + & MPGRS.VPOCHA(NLCEG,INDCEL + 1 )*DYG + & MPGRS.VPOCHA(NLCEG,INDCEL + 2 )*DZG ALCEL = VALCEL + ALCEL * DCEL SCALPA.FRAMG(I1) = ALCEL ENDDO C VALCEL = MPVITP.VPOCHA(NLCEG, 1) ALCEL = MELALV.VPOCHA(NLCEG, 1) DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG + & MPGRV.VPOCHA(NLCEG, 2)*DYG + & MPGRV.VPOCHA(NLCEG, 3)*DZG UXG = VALCEL + ALCEL * DCEL C VALCEL = MPVITP.VPOCHA(NLCEG, 2) ALCEL = MELALV.VPOCHA(NLCEG, 2) DCEL = MPGRV.VPOCHA(NLCEG, 4)*DXG + & MPGRV.VPOCHA(NLCEG, 5)*DYG + & MPGRV.VPOCHA(NLCEG, 6)*DZG UYG = VALCEL + ALCEL * DCEL C VALCEL = MPVITP.VPOCHA(NLCEG, 3) ALCEL = MELALV.VPOCHA(NLCEG, 3) DCEL = MPGRV.VPOCHA(NLCEG, 7)*DXG + & MPGRV.VPOCHA(NLCEG, 8)*DYG + & MPGRV.VPOCHA(NLCEG, 9)*DZG UZG = VALCEL + ALCEL * DCEL C UNG = UXG * CNX + UYG * CNY + UZG * CNZ UTG = UXG * CTX + UYG * CTY + UZG * CTZ UVG = UXG * CVX + UYG * CVY + UZG * CVZ C C********** Si l'on fait pas de prediction, ce n'est pas necessaire de C controller la positivite' de la pression et de la densité; elle C est déjà garantie par la proprieté LED de limiteur; mais, vu C que le limiteur n'est pas calculé ici, mais dans un autre C operateur, on le fait C LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2 C IF(LOGI1)THEN C C************* Premier ordre en espace local C ROG = MPROC.VPOCHA(NLCEG,1) ROD = ROG PG = MPPC.VPOCHA(NLCEG,1) PD = PG UNG = MPVITC.VPOCHA(NLCEG,1)*CNX + & MPVITC.VPOCHA(NLCEG,2)*CNY + & MPVITC.VPOCHA(NLCEG,3)*CNZ UTG = MPVITC.VPOCHA(NLCEG,1)*CTX + & MPVITC.VPOCHA(NLCEG,2)*CTY + & MPVITC.VPOCHA(NLCEG,3)*CTZ UVG = MPVITC.VPOCHA(NLCEG,1)*CVX + & MPVITC.VPOCHA(NLCEG,2)*CVY + & MPVITC.VPOCHA(NLCEG,3)*CVZ UND = -1.0D0 * UNG UTD = UTG UVD = UVG DO I1 = 1, NESP, 1 FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1) FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1) ENDDO DO I1 = 1, NSCA, 1 SCALPA.FRAMG(I1) = MPSCAC.VPOCHA(NLCEG,I1) SCALPA.FRAMD(I1) = SCALPA.FRAMG(I1) ENDDO ELSE C C********** Son etat droite C ROD = ROG PD = PG UND = -1.0D0 * UNG UTD = UTG UVD = UVG DO I1 = 1, NESP, 1 FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1) ENDDO DO I1 = 1, NSCA, 1 SCALPA.FRAMD(I1) = SCALPA.FRAMG(I1) ENDDO ENDIF C C************* Fin cas mur C ELSE C C************* Etat gauche C VALCEL = MPROP.VPOCHA(NLCEG, 1) ALCEL = MELALR.VPOCHA(NLCEG, 1) DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG + & MPGRR.VPOCHA(NLCEG, 2)*DYG + & MPGRR.VPOCHA(NLCEG, 3)*DZG ROG = VALCEL + ALCEL * DCEL C VALCEL = MPPP.VPOCHA(NLCEG, 1) ALCEL = MELALP.VPOCHA(NLCEG, 1) DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG + & MPGRP.VPOCHA(NLCEG, 2)*DYG + & MPGRP.VPOCHA(NLCEG, 3)*DZG PG = VALCEL + ALCEL * DCEL C LOGI2 = .FALSE. SUMY = 0.0D0 DO I1 = 1, NESP, 1 INDCEL = 3 * (I1-1) + 1 VALCEL = MPYP.VPOCHA(NLCEG,I1) ALCEL = MELALY.VPOCHA(NLCEG, I1) DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG + & MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG+ & MPGRY.VPOCHA(NLCEG,INDCEL + 2 )*DZG ALCEL = VALCEL + ALCEL * DCEL SUMY = SUMY + ALCEL LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0) FRAMAS.FRAMG(I1) = ALCEL ENDDO LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0) C DO I1 = 1, NSCA, 1 INDCEL = 3 * (I1-1) + 1 VALCEL = MPSCAP.VPOCHA(NLCEG,I1) ALCEL = MELALS.VPOCHA(NLCEG, I1) DCEL = MPGRS.VPOCHA(NLCEG, INDCEL)*DXG + & MPGRS.VPOCHA(NLCEG,INDCEL + 1 )*DYG+ & MPGRS.VPOCHA(NLCEG,INDCEL + 2 )*DZG ALCEL = VALCEL + ALCEL * DCEL SCALPA.FRAMG(I1) = ALCEL ENDDO C VALCEL = MPVITP.VPOCHA(NLCEG, 1) ALCEL = MELALV.VPOCHA(NLCEG, 1) DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG + & MPGRV.VPOCHA(NLCEG, 2)*DYG + & MPGRV.VPOCHA(NLCEG, 3)*DZG UXG = VALCEL + ALCEL * DCEL C VALCEL = MPVITP.VPOCHA(NLCEG, 2) ALCEL = MELALV.VPOCHA(NLCEG, 2) DCEL = MPGRV.VPOCHA(NLCEG, 4)*DXG + & MPGRV.VPOCHA(NLCEG, 5)*DYG + & MPGRV.VPOCHA(NLCEG, 6)*DZG UYG = VALCEL + ALCEL * DCEL C VALCEL = MPVITP.VPOCHA(NLCEG, 3) ALCEL = MELALV.VPOCHA(NLCEG, 3) DCEL = MPGRV.VPOCHA(NLCEG, 7)*DXG + & MPGRV.VPOCHA(NLCEG, 8)*DYG + & MPGRV.VPOCHA(NLCEG, 9)*DZG UZG = VALCEL + ALCEL * DCEL C UNG = UXG * CNX + UYG * CNY + UZG * CNZ UTG = UXG * CTX + UYG * CTY + UZG * CTZ UVG = UXG * CVX + UYG * CVY + UZG * CVZ C C C********** Positivite C LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2 C IF(LOGI1)THEN C C************* Premier ordre en espace local C ROG = MPROC.VPOCHA(NLCEG,1) PG = MPPC.VPOCHA(NLCEG,1) UNG = MPVITC.VPOCHA(NLCEG,1)*CNX + & MPVITC.VPOCHA(NLCEG,2)*CNY + & MPVITC.VPOCHA(NLCEG,3)*CNZ UTG = MPVITC.VPOCHA(NLCEG,1)*CTX + & MPVITC.VPOCHA(NLCEG,2)*CTY + & MPVITC.VPOCHA(NLCEG,3)*CTZ UVG = MPVITC.VPOCHA(NLCEG,1)*CVX + & MPVITC.VPOCHA(NLCEG,2)*CVY + & MPVITC.VPOCHA(NLCEG,3)*CVZ DO I1 = 1, NESP FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1) ENDDO DO I1 = 1, NSCA, 1 SCALPA.FRAMG(I1) = MPSCAC.VPOCHA(NLCEG,I1) ENDDO ENDIF C C********** Etat droite C VALCEL = MPROP.VPOCHA(NLCED, 1) ALCEL = MELALR.VPOCHA(NLCED, 1) DCEL = MPGRR.VPOCHA(NLCED, 1)*DXD + & MPGRR.VPOCHA(NLCED, 2)*DYD + & MPGRR.VPOCHA(NLCED, 3)*DZD ROD = VALCEL + ALCEL * DCEL C VALCEL = MPPP.VPOCHA(NLCED, 1) ALCEL = MELALP.VPOCHA(NLCED, 1) DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD + & MPGRP.VPOCHA(NLCED, 2)*DYD + & MPGRP.VPOCHA(NLCED, 3)*DZD PD = VALCEL + ALCEL * DCEL C LOGI2 = .FALSE. SUMY = 0.0D0 DO I1 = 1, NESP, 1 INDCEL = 3 * (I1-1) + 1 VALCEL = MPYP.VPOCHA(NLCED,I1) ALCEL = MELALY.VPOCHA(NLCED, I1) DCEL = MPGRY.VPOCHA(NLCED, INDCEL)*DXD + & MPGRY.VPOCHA(NLCED,INDCEL + 1 )*DYD + & MPGRY.VPOCHA(NLCED,INDCEL + 2 )*DZD ALCEL = VALCEL + ALCEL * DCEL SUMY = SUMY + ALCEL LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0) FRAMAS.FRAMD(I1) = ALCEL ENDDO LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0) C DO I1 = 1, NSCA, 1 INDCEL = 3 * (I1-1) + 1 VALCEL = MPSCAP.VPOCHA(NLCED,I1) ALCEL = MELALS.VPOCHA(NLCED, I1) DCEL = MPGRS.VPOCHA(NLCED, INDCEL)*DXD + & MPGRS.VPOCHA(NLCED,INDCEL + 1 )*DYD+ & MPGRS.VPOCHA(NLCED,INDCEL + 2 )*DZD ALCEL = VALCEL + ALCEL * DCEL SCALPA.FRAMD(I1) = ALCEL ENDDO C VALCEL = MPVITP.VPOCHA(NLCED, 1) ALCEL = MELALV.VPOCHA(NLCED, 1) DCEL = MPGRV.VPOCHA(NLCED, 1)*DXD + & MPGRV.VPOCHA(NLCED, 2)*DYD + & MPGRV.VPOCHA(NLCED, 3)*DZD UXD = VALCEL + ALCEL * DCEL C VALCEL = MPVITP.VPOCHA(NLCED, 2) ALCEL = MELALV.VPOCHA(NLCED, 2) DCEL = MPGRV.VPOCHA(NLCED, 4)*DXD + & MPGRV.VPOCHA(NLCED, 5)*DYD + & MPGRV.VPOCHA(NLCED, 6)*DZD UYD = VALCEL + ALCEL * DCEL C VALCEL = MPVITP.VPOCHA(NLCED, 3) ALCEL = MELALV.VPOCHA(NLCED, 3) DCEL = MPGRV.VPOCHA(NLCED, 7)*DXD + & MPGRV.VPOCHA(NLCED, 8)*DYD + & MPGRV.VPOCHA(NLCED, 9)*DZD UZD = VALCEL + ALCEL * DCEL C UND = UXD * CNX + UYD * CNY + UZD * CNZ UTD = UXD * CTX + UYD * CTY + UZD * CTZ UVD = UXD * CVX + UYD * CVY + UZD * CVZ C C********** Positivite C LOGI1 = (PD .LT. 0.0D0) .OR. (ROD .LT. 0.0D0) .OR. LOGI2 C IF(LOGI1)THEN C C************* Premier ordre en espace local C ROD = MPROC.VPOCHA(NLCED,1) PD = MPPC.VPOCHA(NLCED,1) UND = MPVITC.VPOCHA(NLCED,1)*CNX + & MPVITC.VPOCHA(NLCED,2)*CNY + & MPVITC.VPOCHA(NLCED,3)*CNZ UTD = MPVITC.VPOCHA(NLCED,1)*CTX + & MPVITC.VPOCHA(NLCED,2)*CTY + & MPVITC.VPOCHA(NLCED,3)*CTZ UVD = MPVITC.VPOCHA(NLCED,1)*CVX + & MPVITC.VPOCHA(NLCED,2)*CVY + & MPVITC.VPOCHA(NLCED,3)*CVZ DO I1 = 1, NESP FRAMAS.FRAMD(I1) = MPYC.VPOCHA(NLCED,I1) ENDDO DO I1 = 1, NSCA, 1 SCALPA.FRAMD(I1) = MPSCAC.VPOCHA(NLCED,I1) ENDDO ENDIF ENDIF C C******** Les MELVALs C MELRO.VELCHE(1,NLCF) = ROG MELRO.VELCHE(3,NLCF) = ROD MELP.VELCHE(1,NLCF) = PG MELP.VELCHE(3,NLCF) = PD MELVUN.VELCHE(1,NLCF) = UNG MELVUN.VELCHE(3,NLCF) = UND MELVUT.VELCHE(1,NLCF) = UTG MELVUT.VELCHE(3,NLCF) = UTD MELVUV.VELCHE(1,NLCF) = UVG MELVUV.VELCHE(3,NLCF) = UVD MELVNX.VELCHE(1,NLCF) = CNX MELVNY.VELCHE(1,NLCF) = CNY MELVNZ.VELCHE(1,NLCF) = CNZ MELT1X.VELCHE(1,NLCF) = CTX MELT1Y.VELCHE(1,NLCF) = CTY MELT1Z.VELCHE(1,NLCF) = CTZ MELT2X.VELCHE(1,NLCF) = CVX MELT2Y.VELCHE(1,NLCF) = CVY MELT2Z.VELCHE(1,NLCF) = CVZ DO I1 = 1, NESP, 1 MELVA1 = MCHAMY.IELVAL(I1) MELVA1.VELCHE(1,NLCF) = FRAMAS.FRAMG(I1) MELVA1.VELCHE(3,NLCF) = FRAMAS.FRAMD(I1) ENDDO DO I1 = 1, NSCA, 1 MELVA1 = MCHAMS.IELVAL(I1) MELVA1.VELCHE(1,NLCF) = SCALPA.FRAMG(I1) MELVA1.VELCHE(3,NLCF) = SCALPA.FRAMD(I1) ENDDO ENDDO C C**** Desactivation des SEGMENTs C SEGDES IPT1 SEGDES IPT2 C C**** Le MPOVALs 'Prediction' sont detruits (si existentes) C IF(LOGTEM)THEN SEGSUP MPROP SEGSUP MPVITP SEGSUP MPPP IF(LYC) SEGSUP MPYP IF(LSCAC) SEGSUP MPSCAP SEGDES MPGAMC ENDIF C SEGDES MPROC SEGDES MPGRR SEGDES MELALR SEGDES MPVITC SEGDES MPGRV SEGDES MELALV SEGDES MPPC SEGDES MPGRP SEGDES MELALP IF(LYC)THEN SEGDES MPYC SEGDES MPGRY SEGDES MELALY SEGDES MPYC DO I1 = 1, NESP, 1 MELVA1 = MCHAMY.IELVAL(I1) SEGDES MELVA1 ENDDO SEGDES MCHAMY SEGSUP FRAMAS ENDIF IF(LSCAC)THEN SEGDES MPSCAC SEGDES MPGRS SEGDES MELALS DO I1 = 1, NSCA, 1 MELVA1 = MCHAMS.IELVAL(I1) SEGDES MELVA1 ENDDO SEGDES MCHAMS SEGSUP SCALPA ENDIF SEGDES MPNORM C SEGDES MELRO SEGDES MELP SEGDES MELVUN SEGDES MELVUT SEGDES MELVUV SEGDES MELVNX SEGDES MELVNY SEGDES MELVNZ SEGDES MELT1X SEGDES MELT1Y SEGDES MELT1Z SEGDES MELT2X SEGDES MELT2Y SEGDES MELT2Z C C**** Destruction du MELNTI correspondance local/global C SEGSUP MLENT1 C 9999 CONTINUE C RETURN END
© Cast3M 2003 - Tous droits réservés.
Mentions légales