Sources Esope | Cast3M

konfl2
C KONFL2    SOURCE    OF166741  24/12/13    21:16:24     12097          
      SUBROUTINE KONFL2(LOGME,LOGSCA,INDMET,NORDP1,
     &          IROF,IVITF,IPF,IFRMAF,ISCAF,PROPHY,
     &          ICHPSU,ICHPDI,
     &          MELEMC,MELEMF,MELEFE,MELLIM,
     &          ICHFLU,DT,
     &          LOGNC,LOGAN,MESERR)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CKON3
C
C DESCRIPTION       :  Voir CKON
C
C                      Cas deux dimensions, gaz "thermally perfect"
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
C
C************************************************************************
C
C
C APPELES (Outils
C          CASTEM)  :  KRIPAD, LICHT
C
C APPELES (Calcul)  :  FVLHTE, FCOLTE
C
C
C************************************************************************
C
C ENTREES
C
C
C  1) PARAMETRES
C
C     LOGME   : (LOGICAL);  .TRUE.  -> MULTI-ESPECES
C                           .FALSE. -> MONO-ESPECE
C
C     LOGSCA  : (LOGICAL);  .TRUE.  -> scalaires passifs
C                           .FALSE. -> no scalaires passifs
C
C     INDMET  :     1        VLH (van Leer Hanel FVS)
C
C                   2        Colella-Glaz FDS
C
C      NORDP1 : (ordre des polynoms cv(T)) + 1
C
C  2) Pointeurs des MCHAMLs
C
C     IROF    : MCHAML sur "FACEL"  contenant la masse volumique
C     ("gauche" et "droite");
C
C     IVITF   : MCHAML sur "FACEL"  contenant la vitesse dans le repaire
C               local (n,t) et les cosinus directeurs des repaire local;
C
C     IPF     : MCHAML sur "FACEL"  contenant la pression;
C
C     IFRAMAF : MCHAML sur "FACEL",  contenant les fractions massiques
C               si LOGME = .TRUE.;
C               LOGME = .FALSE. -> IFRAMAF = 0
C
C     ISCAF   : MCHAML sur "FACEL",  contenant les scalires passifs a
C               transporter (ou ISCAF = 0)
C
C  3) Pointeur sur la table des proprietes de gaz
C
C     PROPHY
C
C  4) Pointeurs de CHPOINTs de la table DOMAINE
C
C     ICHPSU  : CHPOINT "FACE" contenant la surface des faces
C
C     ICHPDI  : CHPOINT "CENTRE" contenant le diametre minimum
C               de chaque element
C
C
C  5) Pointeurs de MELEME de la table DOMAINE
C
C     MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
C
C     MELEMF  : MELEME 'FACE' du SPG des FACES
C
C     MELEFE  : MELEME 'FACEL' du connectivité FACES -> ELEM
C
C     MELLIM  : MAILLAGE where the flux (or residu) will not be found
C
C SORTIES (il faudrait dire E/S)
C
C     ICHFLU     : pointeurs de CHPOINTs "FACE" des flux aux interfaces:
C
C     DT      : pas de temps pour le respect de la CFL-like condition
C               DT &lt; DIAMEL /2 /max(Lambda_i)
C               En maillage regulier cette condition garantie la
C               non-interaction des ondes
C
C
C     LOGNC   : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
C               dans pour la solution du probleme Riemann n'a pas bien
C               marchéee
C
C     LOGAN   : (LOGICAL): si .TRUE. une anomalie à été detectée
C
C     MESERR  : pour l'ecriture des messages d'erreurs
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE : 22.12.98 Creation
C
C              08.02.00 Ajout du flux numerique de Colella-Glaz
C
C              21.02.00 Ajout de transport de scalaires passifs
C
C              03.12.03 Ajout le MAILLIM comme le l'entree
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         GAMMA \in (1,3)
C                                         Y \in (0,1)
C       Si non il faut le faire!!!
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      INTEGER I1, I2
     &       ,INDMET, NORDP1
     &       ,IROF,IVITF,IPF,ISCAF,IFRMAF
     &       ,ICHPSU,ICHPDI,MELEMC,MELEMF,MELEFE
     &       ,IGEOMC,IGEOMF
     &       ,ICHFLU
     &       ,NESP, NFAC, NSCA, NMA, NMA2
     &       ,NLCF, NGCEG,  NGCED, NLCEG, NLCED
     &       ,NGCF, NLCF1, SPG1, SPG2, NLFL
      REAL*8 DT, UNSDT, CELLT
     &     , ROG, UNG, UTG, PG, GAMG, TG
     &     , ROD, UND, UTD, PD, GAMD, TD
     &     , SURF,CNX, CNY, CTX , CTY
     &     , CELL, DIAMG, DIAMD, DIAM
     &     , ASON, LAMBDA
     &     , Y1G, Y1D, RTOTG, RTOTD, CVTOTG, CVTOTD, ETHERG, ETHERD
     &     , YG, YD, FUNTG, FUNTD, DCV, XST
      LOGICAL LOGME, LOGSCA, LOGNC, LOGAN
      CHARACTER*(40) MESERR
      CHARACTER*(8) TYPE
C
C**** Segment des proprietes du gaz
C
      SEGMENT PROPHY
         REAL*8 ACV(NORDP1,NESP+1), R(NESP+1), H0K(NESP+1)
     &         ,ACVTOG(NORDP1), ACVTOD(NORDP1)
      ENDSEGMENT
C
C**** Les Includes
C
C
C**** LES INCLUDES
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMCHAML
      POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
     &         MELT1X.MELVAL, MELT1Y.MELVAL
      POINTEUR MELVUN.MELVAL, MELVUT.MELVAL
      POINTEUR MELRO.MELVAL, MELP.MELVAL
-INC SMCHPOI
      POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
     &        , MPOFLU.MPOVAL
      POINTEUR MCHAMY.MCHAML,  MCHAMS.MCHAML
-INC SMELEME
      POINTEUR MELLIM.MELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MLELIM.MLENTI
C
C**** Les fractiones massiques.
C
      SEGMENT FRAMAS
         REAL*8 YET(NMA)
      ENDSEGMENT
      POINTEUR FRAMAG.FRAMAS, FRAMAD.FRAMAS, SCAG.FRAMAS, SCAD.FRAMAS
C
C**** Les flux aux interface dans le repaire (n,t)
C
      SEGMENT IFLUX
         REAL*8 FLUX(NMA2)
      ENDSEGMENT
      POINTEUR IFLU1.IFLUX, IFLU2.IFLUX
C
C**** KRIPAD pour la correspondance global/local des conditions limits
C
      CALL KRIPAD(MELLIM,MLELIM)
C
C**** Initialisation des MCHAMLs
C
C**** Masse volumique
C
      MCHEL1 = IROF
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELRO = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C**** Pression
C
      MCHEL1 = IPF
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELP = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C**** Vitesse et cosinus directeurs du repere (n,t)
C
      MCHEL1 = IVITF
      SEGACT MCHEL1
C
C**** La vitesse a comme SPG MELEFE
C     Le cosinus directeurs ont comme SPG MELEMF
C
C     MCHAM1 -> Cosinus directeurs
C     MCHAM2 -> Vitesse
C
      SPG1 = MCHEL1.IMACHE(1)
      SPG2 = MCHEL1.IMACHE(2)
      IF((SPG1 .EQ. MELEMF)  .AND. (SPG2 .EQ. MELEFE))THEN
         MCHAM1 = MCHEL1.ICHAML(1)
         MCHAM2 = MCHEL1.ICHAML(2)
      ELSEIF((SPG1 .EQ. MELEFE)  .AND. (SPG2 .EQ. MELEMF))THEN
         MCHAM1 = MCHEL1.ICHAML(2)
         MCHAM2 = MCHEL1.ICHAML(1)
      ELSE
         LOGAN = .TRUE.
         GOTO 9999
      ENDIF
      SEGACT MCHAM1
      MELVNX = MCHAM1.IELVAL(1)
      MELVNY = MCHAM1.IELVAL(2)
      MELT1X = MCHAM1.IELVAL(3)
      MELT1Y = MCHAM1.IELVAL(4)
      SEGDES MCHAM1
      SEGACT MCHAM2
      MELVUN = MCHAM2.IELVAL(1)
      MELVUT = MCHAM2.IELVAL(2)
      SEGDES MCHAM2
      SEGDES MCHEL1
C
C**** Fractions massiques
C
      IF(LOGME)THEN
         MCHEL1 = IFRMAF
         SEGACT MCHEL1
         MCHAMY = MCHEL1.ICHAML(1)
         SEGACT MCHAMY
C
C******* Numero d'especes dans les equations d'Euler
C
         NESP =  MCHAMY.IELVAL(/1)
         DO I1 = 1, NESP
            MELVA1 = MCHAMY.IELVAL(I1)
            SEGACT MELVA1
         ENDDO
         NMA = NESP
         SEGINI FRAMAG
         SEGINI FRAMAD
         SEGDES MCHEL1
      ELSE
C
C******* Definition minimale de YET, necessaire pour transmetre YET aux
C        subroutines FORTRAN qui calculent les flux
C
         NMA = 1
         SEGINI FRAMAG
         SEGINI FRAMAD
         NESP = 0
      ENDIF
C
C**** Scalaires passifs
C
      IF(LOGSCA)THEN
         MCHEL1 = ISCAF
         SEGACT MCHEL1
         MCHAMS = MCHEL1.ICHAML(1)
         SEGACT MCHAMS
         NSCA =  MCHAMS.IELVAL(/1)
         DO I1 = 1, NSCA, 1
            MELVA1 = MCHAMS.IELVAL(I1)
            SEGACT MELVA1
         ENDDO
         NMA = NSCA
         SEGINI SCAG
         SEGINI SCAD
         SEGDES MCHEL1
      ELSE
C
C******* Definition minimale de YET, necessaire pour transmetre YET aux
C        subroutines FORTRAN qui calculent les flux
C
         NMA = 1
         SEGINI SCAG
         SEGINI SCAD
         NSCA= 0
      ENDIF
C
C**** Initialisation des MELEMEs
C
C     'CENTRE', 'FACEL'
C
      IPT2 = MELEFE
      SEGACT IPT2
      NFAC = IPT2.NUM(/2)
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(MELEMC,MLENT1)
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
C**** KRIPAD pour la correspondance global/local de 'FACE'
C
      CALL KRIPAD(MELEMF,MLENT2)
C
C**** Initialisation de flux
C
      NMA2 = 4 + NESP + NSCA
      SEGINI IFLU1
      SEGINI IFLU2
C
C**** IFLU2 = segment de travail en FLUVLH; c'est plus rapide le definir ici
C
C
C**** CHPOINTs de la table DOMAINE
C
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
C
C**** LICHT active les MPOVALs en *MOD
C
C     i.e.
C
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVDI*MOD
C
C
C**** Les FLUX aux face
C
C     La densité
C
      CALL LICHT(ICHFLU,MPOFLU,TYPE,IGEOMF)
C
C     SEGACT MPOFLU*MOD
C
C**** Activation des MCHAMLs
C
      SEGACT MELRO
      SEGACT MELP
      SEGACT MELVUN
      SEGACT MELVUT
      SEGACT MELVNX
      SEGACT MELVNY
      SEGACT MELT1X
      SEGACT MELT1Y
C
C**** Initialisation de 1/DT
C
      UNSDT = 0.0D0
C
C**** BOUCLE SUR FACEL pour le calcul du FLUX
C
      DO NLCF = 1, NFAC
C
C******* NLCF  = numero local du centre de facel
C        NGCF  = numero global du centre de facel
C        NLCF1 = numero local 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 = IPT2.NUM(1,NLCF)
         NGCED = IPT2.NUM(3,NLCF)
         NGCF  = IPT2.NUM(2,NLCF)
         NLCF1 = MLENT2.LECT(NGCF)
         NLCEG = MLENT1.LECT(NGCEG)
         NLCED = MLENT1.LECT(NGCED)
         NLFL  = MLELIM.LECT(NGCF)
C
C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
C
         IF(NLCF .NE. NLCF1)THEN
            MESERR = 'Il ne faut pas jouer avec la console.   '
            LOGAN = .TRUE.
            GOTO 9999
         ENDIF
         IF(NLFL .EQ. 0)THEN
C
C******* Recuperation des Etats "gauche" et "droite"
C
         ROG = MELRO.VELCHE(1,NLCF)
         UNG = MELVUN.VELCHE(1,NLCF)
         UTG = MELVUT.VELCHE(1,NLCF)
         PG = MELP.VELCHE(1,NLCF)
C
         ROD = MELRO.VELCHE(3,NLCF)
         UND = MELVUN.VELCHE(3,NLCF)
         UTD = MELVUT.VELCHE(3,NLCF)
         PD = MELP.VELCHE(3,NLCF)
C
         CNX = MELVNX.VELCHE(1,NLCF)
         CNY = MELVNY.VELCHE(1,NLCF)
         CTX = MELT1X.VELCHE(1,NLCF)
         CTY = MELT1Y.VELCHE(1,NLCF)
C
C******* Le fractiones massiques, R et ACVTO
C
         RTOTG = 0.0D0
         RTOTD = 0.0D0
         CVTOTG = 0.0D0
         CVTOTD = 0.0D0
         Y1G = 1.0D0
         Y1D = 1.0D0
         ETHERG = 0.0D0
         ETHERD = 0.0D0
         DO I1= 1, NORDP1, 1
            PROPHY.ACVTOG(I1) =0.0D0
            PROPHY.ACVTOD(I1) =0.0D0
         ENDDO
         DO I1 = 1, NESP, 1
            MELVA1 = MCHAMY.IELVAL(I1)
            YG = MELVA1.VELCHE(1,NLCF)
            YD = MELVA1.VELCHE(3,NLCF)
            Y1G = Y1G - YG
            Y1D = Y1D - YD
            FRAMAG.YET(I1) = YG
            FRAMAD.YET(I1) = YD
            RTOTG = RTOTG + YG * PROPHY.R(I1)
            RTOTD = RTOTD + YD * PROPHY.R(I1)
            DO I2 = 1, NORDP1
               PROPHY.ACVTOG(I2) = PROPHY.ACVTOG(I2) +
     &                             YG * PROPHY.ACV(I2,I1)
               PROPHY.ACVTOD(I2) = PROPHY.ACVTOD(I2) +
     &                             YD * PROPHY.ACV(I2,I1)
            ENDDO
         ENDDO
         RTOTG = RTOTG + Y1G * PROPHY.R(NESP+1)
         RTOTD = RTOTD + Y1D * PROPHY.R(NESP+1)
         DO I2 = 1, NORDP1, 1
            PROPHY.ACVTOG(I2) = PROPHY.ACVTOG(I2) +
     &         Y1G * PROPHY.ACV(I2,NESP+1)
            PROPHY.ACVTOD(I2) = PROPHY.ACVTOD(I2) +
     &         Y1D * PROPHY.ACV(I2,NESP+1)
         ENDDO
         TG = PG / (ROG * RTOTG)
         TD = PD / (ROD * RTOTD)
         FUNTG = 1.0D0
         FUNTD = 1.0D0
         CVTOTG = PROPHY.ACVTOG(1)
         ETHERG = CVTOTG * TG
         CVTOTD = PROPHY.ACVTOD(1)
         ETHERD = CVTOTD * TD
         DO I1 = 2, NORDP1, 1
            FUNTG = FUNTG * TG
            FUNTD = FUNTD * TD
            DCV = PROPHY.ACVTOG(I1) * FUNTG
            CVTOTG = CVTOTG + DCV
            ETHERG = ETHERG + DCV * TG / I1
            DCV = PROPHY.ACVTOD(I1) * FUNTD
            CVTOTD = CVTOTD + DCV
            ETHERD = ETHERD + DCV * TD / I1
         ENDDO
         GAMG = (CVTOTG + RTOTG) /CVTOTG
         GAMD = (CVTOTD + RTOTD) /CVTOTD
C
C******* Les scalaires passifs
C
         DO I1 = 1, NSCA, 1
            MELVA1 = MCHAMS.IELVAL(I1)
            SCAG.YET(I1) = MELVA1.VELCHE(1,NLCF)
            SCAD.YET(I1) = MELVA1.VELCHE(3,NLCF)
         ENDDO
C
C******* On a defini (ROg,ROUNg,ROUTg,Pg,(Yg)), (ROd,ROUNd,ROUTd,Pd,(Yd))
C        et on a déjà verifié ROg, ROd, Pg, Pd > 0 et 0&lt;Y_i&lt;1
C
C
C******* Calcul du flux aux interfaces
C
         IF(INDMET .EQ. 1)THEN
C
C********** VLH FVS
C
            CALL FVLHTE(NESP,NSCA,
     &                  GAMG,ROG,PG,UNG,UTG,ETHERG,
     &                  GAMD,ROD,PD,UND,UTD,ETHERD,
     &                  FRAMAG.YET,FRAMAD.YET,
     &                  SCAG.YET,SCAD.YET,
     &                  IFLU1.FLUX,IFLU2.FLUX,
     &                  CELLT)
         ELSEIF(INDMET .EQ. 2)THEN
C
C******* Colella-Glaz FDS (avec Entropy-fix)
C
            XST=0.0D0
            CALL FCOLTE(NESP,NSCA,NORDP1,XST,
     &               GAMG,RTOTG,PROPHY.ACVTOG,ROG,PG,TG,UNG,UTG,ETHERG,
     &               GAMD,RTOTD,PROPHY.ACVTOD,ROD,PD,TD,UND,UTD,ETHERD,
     &               FRAMAG.YET,FRAMAD.YET,
     &               SCAG.YET,SCAD.YET,
     &               IFLU1.FLUX,CELLT,
     &               LOGNC,MESERR,LOGAN)
         ENDIF
C
         IF(LOGAN) GOTO 9999
         IF(LOGNC) GOTO 9999
C
C******* Ecriture des flux
C
C        FLUX(1) = RO Un              RO Un
C        FLUX(2) = RO Un Un + P   ->  RO Un Ux + P CNX
C        FLUX(3) = RO Un Ut       ->  RO Un Uy + P CNY
C        FLUX(4) = RO Un Et           RO Un Et
C
         SURF   = MPOVSU.VPOCHA(NLCF,1)
         MPOFLU.VPOCHA(NLCF,1) =
     &        (IFLU1.FLUX(1) * SURF )
         MPOFLU.VPOCHA(NLCF,2) =
     &        ((IFLU1.FLUX(2)*CNX+IFLU1.FLUX(3)*CTX) * SURF)
         MPOFLU.VPOCHA(NLCF,3) =
     &        ((IFLU1.FLUX(2)*CNY+IFLU1.FLUX(3)*CTY) * SURF)
         MPOFLU.VPOCHA(NLCF,4) =
     &        (IFLU1.FLUX(4) * SURF)
         DO I1 = 1, NESP, 1
            MPOFLU.VPOCHA(NLCF,4+I1)=IFLU1.FLUX(4+I1)
     &           * SURF
         ENDDO
         DO I1 = 1, NSCA, 1
            MPOFLU.VPOCHA(NLCF,4+NESP+I1)=IFLU1.FLUX(4+I1+NESP)
     &           * SURF
         ENDDO
C
C******* Calcul du pas du temps (CFL)
C
C******  a) etat a l'interface
C
         DIAMG = MPOVDI.VPOCHA(NLCEG,1)
         DIAMD = MPOVDI.VPOCHA(NLCED,1)
         DIAM = MIN(DIAMG,DIAMD)
         CELL  = 1.0D0/DIAM/CELLT
         IF(CELL .GT. UNSDT)THEN
            UNSDT  =  CELL
         ENDIF
C
C******  b) etat gauche
C
         ASON    = SQRT(GAMG*PG/ROG)
         LAMBDA = ABS(UNG) + ASON
         CELL = LAMBDA / DIAM
         IF(CELL .GT. UNSDT)THEN
            UNSDT  =  CELL
         ENDIF
C
C******  C) etat droite
C
         ASON    = SQRT(GAMD*PD/ROD)
         LAMBDA = ABS(UND) + ASON
         CELL = LAMBDA / DIAM
         IF(CELL .GT. UNSDT)THEN
            UNSDT  =  CELL
         ENDIF
        ENDIF
C
C**** Fin boucle sur FACEL
C
      ENDDO
C
C**** Pas du temps (condition de non interaction en 1D)
C
      DT = 0.5D0 / UNSDT
C
C**** Desactivation des segments et
C     on detruit les MCHAMLs
C
C
C**** SEGSUP FRAMAG
C     SEGSUP FRAMAD
C
C     meme si LOGME = .FALSE.
C
      SEGDES FRAMAG
      SEGDES FRAMAD
C
      SEGDES MLENT1
      SEGDES MLENT2
      SEGDES IPT2
C
      SEGDES IFLU1
      SEGDES IFLU2
C
      SEGDES MPOVSU
      SEGDES MPOVDI
C
      SEGDES MPOFLU
C
      SEGDES MELRO
      SEGDES MELP
      SEGDES MELVUN
      SEGDES MELVUT
      SEGDES MELVNX
      SEGDES MELVNY
      SEGDES MELT1X
      SEGDES MELT1Y
c      SEGDES MELLIM
      SEGDES MLELIM
C
      IF(LOGME) THEN
         DO I1 = 1, NESP
            MELVA1 = MCHAMY.IELVAL(I1)
            SEGDES MELVA1
         ENDDO
C
         SEGDES MCHAMY
      ENDIF
      IF(LOGSCA) THEN
         DO I1 = 1, NSCA, 1
            MELVA1 = MCHAMS.IELVAL(I1)
            SEGDES MELVA1
         ENDDO
C
         SEGDES MCHAMS
      ENDIF
C
 9999 CONTINUE
C
      RETURN
      END
C