Sources Esope | Cast3M

kon11
C KON11     SOURCE    CB215821  20/11/25    13:31:48     10792          
      SUBROUTINE KON11(IFLU,LOGME)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KON11
C
C DESCRIPTION       :  Subroutine appellée par KONV1
C
C                      Modelisation 2D/3D des equations d'Euler
C
C                      Calcul du flux/residu, DELTAT
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C APPELES (Calcul)  :  KONFL1 (2D, gaz "calorically perfect")
C                      KONFL3 (3D, gaz "calorically perfect")
C
C************************************************************************
C
C ENTRÉES : IFLU    : si EGAL à 1 on veut calculer le flux
C           IFLU    : si EGAL à 2 on veut calculer le residu
C
C           LOGME   : .TRUE./.FALSE. Multi-espece/mono-espece
C
C*** SYNTAXE
C
C  Discrétisation en VF "cell-centered" des équations d'Euler pour
C  un gaz parfait polytropique
C  Inconnues: densités, quantité de mouvement, énergie totale par
C  unité de volumes (variables conservatives)
C
C  RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMONO' MOT1 MOT2 LMOT1
C                  MOD1 MCHAM1 MCHAM2 MCHAM3 MCHAM4 (MAIL1) ;
C
C  or (Bas Mach mono-espece)
C
C  RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMONO' MOT1 MOT2 LMOT1
C                   MOD1 MCHAM1 MCHAM2 MCHAM3 MCHAM4
C                   CHPO1 CHPO2 (MAIL1) ;
C
C  or (multiespece)
C
C  RCHPO1 RFLOT1 = 'KONV' 'VF' 'PERFMULT' MOT1 MOT2 LMOT1
C                   MOD1 MCHAM1 MCHAM2 MCHAM3 MCHAM4 MCHAM5 TABP ;
C
C  ENTREES
C
C  MOT1   : objet de type MOT
C           Il vaut 'RESI' si on veut calculer le résidu
C           Il vaut 'FLUX' si on veut calculer le flux
C
C  MOT2   : objet de type MOT
C           Il indique la méthode de décentrement:
C           'GODUNOV'  = solveur exacte
C           'VANLEER'  = solveur de van Leer
C           'VLH'      = solveur de van Leer Hanel
C           'HUSVL'    = HUS (van Leer + Osher)
C           'HUSVLH'   = HUS (van Leer Hanel + Osher)
C           'AUSMPLUS' = AUSM+
C           'ROE'      = solveur de Roe
C           'SS'       = solveur choc-choc
C           'RUSANOV'  = solveur de Rusanov
C           'AUSMPLM'  = AUSM+ low Mach
C           'RUSANOLM' = solveur de Rusanov low Mach
C           'CENTERED' = solveur centré
C           'ROELM'    = solveur de Roe-Turkel
C           'HLLC'     = solveur 'HLLC'
C           'HLLCLM'   = solveur 'HLLC'-Turkel low Mach
C           'AUSMPUP'  = AUSM+up low Mach
C
C
C LMOT1   : objet de type LISTMOTS
C           Noms de composantes du résultat (RCHPO1)
C           Il contient dans l'ordre suivant: le noms de la densité,
C           du momentum, de l'énergie totale par unité de volume
C           (le densité des constituants qui sont dans TABP.'ESPEULE'
C            dans le cas multi-espece)
C
C MOD1    : objet modele de type Navier_Stokes
C
C MCHAM1  : MCHAML contenant la masse volumique, qui a comme
C           SPG (support géométrique) l'indice 'FACEL' de
C           MOD1 (une composante, 'SCAL')
C
C MCHAM2  : MCHAML  contenant la vitesse  et les cosinus
C           directeurs du repère locale (n,t) dans le repère
C           global (x,y) (dans le cas 2D 6 composantes:
C           * 'UN' = vitesse normale  (SPG =('DOMA' MOD1 'FACEL'))
C           * 'UT' = vitesse tangentielle  (SPG =('DOMA' MOD1 'FACEL'))
C           * 'NX' = n.x (SPG = 'FACE')
C           * 'NY' = n.y (SPG = 'FACE')
C           * 'TX' = t.x (SPG = 'FACE')
C           * 'TY' = t.y (SPG = 'FACE')).
C
C MCHAM3  : MCHAML (SPG =('DOMA' MOD1 'FACEL')) contenant la pression du
C           gaz (une  seule composante, 'SCAL').
C
C MCHAM4  : MCHAML (SPG =('DOMA' MOD1 'FACEL')) contenant le "gamma" du
C           gaz (une seule composante, 'SCAL').
C
C MCHAM5  : MCHAML (SPG =('DOMA' MOD1 'FACEL')) contenant les fractions
C           massiques du gaz dans le cas multi-espece
C           (noms de composantes dans TABP . 'ESPEULE')
C
C TABP    : la table des proprietés du gaz (voir operateur 'PRIM',
C           modele 'PERFMULT')
C
C CHPO1   : vitesse de cut-off (SPG =('DOMA' MOD1 'CENTRE'), une seule
C           composante, 'SCAL')
C           A donner que pour les solveurs Bas Mach
C
C CHPO2   : deuxieme vitesse de cut-off  (SPG = ('DOMA' MOD1 'CENTRE'),
C           une seule composante, 'SCAL')
C           A donner que pour les solveurs Bas Mach
C
C (MAIL1) : maillage de 'POI1' qui contient de point 'FACE' ou on ne
C           veut pas calculer les flux
C
C SORTIES
C
C
C RCHPO1  : objet de type CHPOINT (composantes =  LMOT1)
C           Residu  si MOT2 = 'RESI' (SPG = ('DOMA' MOD1 'CENTRE'))
C           Flux si MOT2 = 'FLUX'    (SPG = ('DOMA' MOD1 'FACE'))
C
C RFLOT1  : objet de type FLOTTANT
C           Il est le temps caracteristique associé à l'onde la plus
C           rapide.
C
C Remarque
C --------
C
C RCHPO1 est égal à:
C * la derivé temporelle des inconnues si l'option 'RESI' est utilisée
C * la projection du flux convectif sur (DOMA MOD1 'XXNORMAF') si
C   l'option 'FLUX' est utilisée
C
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C              01/05/2005 par T. KLOCZKO (DEN/DM2S/SFME/LTMF)
C              Ajout d'une méthode de décentrement pour les écoulements
C              bas-Mach.
C
C  MOT2 :  'ROELM'  =  schéma de ROE avec préconditionnement Turkel
C
C************************************************************************
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC SMLMOTS
-INC SMCHPOI
-INC SMELEME
      POINTEUR MLMVIT.MLMOTS,MLMOTY.MLMOTS,MLMODT.MLMOTS
C
      INTEGER IFLU, NBMET, INDMET, IRET, INDIC, NBCOMP
     &     ,IDOMA, MELEMC, MELEMF, MELEFE, MELLIM, ICHPSU, ICHPDI
     &     ,ICHPVO, INORM, JGN, JGM
     &     , NC, NESP, ICACCA, IPGAZ
     &     ,ILIINC
     &     ,IROF, IVITF, IPF, IGAMF, IFRMAF, IUINF, IUPRI, IFLIM
     &     ,ICHFLU, ICHRES, INEFMD, ICOND, MMODEL
      PARAMETER (NBMET=16)
      REAL*8 DT
      CHARACTER*8 LMETO(NBMET), TYPE
      CHARACTER*4 MOT
      CHARACTER*(40) MESERR
      LOGICAL LOGNC, LOGAN, LOGME
C
      DATA LMETO/'GODUNOV ','VANLEER','VLH     ',
     &           'HUSVL   ','HUSVLH ','AUSMPLUS','ROE     ',
     &           'SS      ','AUSMPLM','RUSANOV', 'RUSANOLM',
     &           'CENTERED','ROELM  ','HLLC   ','HLLCLM ',
     &           'AUSMPUP ' /
C
C**** Initialisation des variables pour la gestion des erreurs.
C
      LOGNC  = .FALSE.
      LOGAN  = .FALSE.
      MESERR  = '                                        '
C
C**** Metode utilisée
C
      CALL LIRMOT(LMETO,NBMET,INDMET,1)
      IF(IERR .NE. 0)GOTO 9999
      IF(INDMET .EQ. 0)THEN
C
C******** Message d'erreur standard
C         251 2
C         Tentative d'utilisation d'une option non implémentée
C
         CALL ERREUR(251)
      ENDIF
C
C**********************************
C**** Lecture de l'objet MODELE ***
C**********************************
C
      ICOND = 1
      CALL QUETYP(TYPE,ICOND,IRET)
C
      IF(IRET.EQ.0.AND.TYPE.NE.'MMODEL')THEN
        WRITE(IOIMP,*)' On attend un objet MMODEL'
        GOTO 9999
      ENDIF
      CALL LIROBJ('MMODEL',MMODEL,ICOND,IRET)
      IF(IERR.NE.0)GOTO 9999
      CALL LEKMOD(MMODEL,IDOMA,INEFMD)
C     INEFMD inutilisé
      IF(IERR.NE.0)GOTO 9999
C
C**** Centre, FACE, FACEL
C
      CALL LEKTAB(IDOMA,'CENTRE',MELEMC)
      IF(IERR .NE. 0) GOTO 9999
C
      CALL LEKTAB(IDOMA,'FACE',MELEMF)
      IF(IERR .NE. 0) GOTO 9999
C
      CALL LEKTAB(IDOMA,'FACEL',MELEFE)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Lecture du CHPOINT contenant les surfaces des faces.
C
      CALL LEKTAB(IDOMA,'XXSURFAC',ICHPSU)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(ICHPSU, MELEMF, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Lecture du CHPOINT contenant les diametres minimums.
C
      CALL LEKTAB(IDOMA,'XXDIEMIN',ICHPDI)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(ICHPDI, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Lecture du CHPOINT contenant les volumes
C
      CALL LEKTAB(IDOMA,'XXVOLUM',ICHPVO)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(ICHPVO, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Les normales aux faces
C
      IF(IDIM .EQ. 2)THEN
C        Que les normales
         CALL LEKTAB(IDOMA,'XXNORMAF',INORM)
         IF(IERR .NE. 0) GOTO 9999
         JGN = 4
         JGM = 2
         SEGINI MLMVIT
         MLMVIT.MOTS(1) = 'UX  '
         MLMVIT.MOTS(2) = 'UY  '
         CALL QUEPO1(INORM, MELEMF, MLMVIT)
         SEGSUP MLMVIT
         IF(IERR .NE. 0) GOTO 9999
      ELSE
C        Les normales et les tangentes
         TYPE = ' '
         CALL ACMO(IDOMA,'MATROT',TYPE,INORM)
         IF (TYPE .NE. 'CHPOINT ') THEN
            CALL MATRAN(IDOMA,INORM)
            IF(IERR .NE. 0) GOTO 9999
         ENDIF
         JGN = 4
         JGM = 9
         SEGINI MLMVIT
         MLMVIT.MOTS(1) = 'UX  '
         MLMVIT.MOTS(2) = 'UY  '
         MLMVIT.MOTS(3) = 'UZ  '
         MLMVIT.MOTS(4) = 'RX  '
         MLMVIT.MOTS(5) = 'RY  '
         MLMVIT.MOTS(6) = 'RZ  '
         MLMVIT.MOTS(7) = 'MX  '
         MLMVIT.MOTS(8) = 'MY  '
         MLMVIT.MOTS(9) = 'MZ  '
         CALL QUEPO1(INORM, MELEMF, MLMVIT)
         SEGSUP MLMVIT
      ENDIF
C
C********************************
C**** Fin table domaine *********
C********************************
C
C**** On va lire les pointeurs des MCHAMLs
C     Lecture du MCHAML 'FACEL' densité
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IROF,1,IRET)
      IF(IERR.NE.0) GOTO 9999
C
C**** Lecture du MCHAML 'FACEL' vitesse
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IVITF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Lecture du MCHAML 'FACEL' contenant la pression
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IPF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Lecture du MCHAML 'FACEL' contenant les gamma
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,IGAMF,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Si LOGME -> MULTIESPECES
C
      IF(LOGME)THEN
C
C******** Lecture du MCHAML 'FACEL' contenant les fractiones massiques
C
         TYPE='MCHAML  '
         CALL  LIROBJ(TYPE,IFRMAF,1,IRET)
         IF(IERR .NE. 0) GOTO 9999
C
C********** Lecture de la table qui contient le proprieté du gaz
C
         TYPE='TABLE   '
         CALL  LIROBJ(TYPE,IPGAZ,1,IRET)
         IF(IERR .NE. 0) GOTO 9999
         TYPE='       '
         CALL ACMO(IPGAZ,'ESPEULE ',TYPE,MLMOTY)
         IF(IERR .NE. 0) GOTO 9999
         IF(TYPE .NE. 'LISTMOTS')THEN
C
C********** Message d'erreur standard
C           -301 0 %m1:40
C
            MOTERR(1:40) = 'TABGAS = ???                            '
            WRITE(IOIMP,*) MOTERR
C
C********** Message d'erreur standard
C           21 2
C           Données incompatibles
C
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
C
         SEGACT MLMOTY
         NESP = MLMOTY.MOTS(/2)
C
      ELSE
         IFRMAF = 0
         NESP=0
      ENDIF
C
C**** La list des inconnues
C
      TYPE='LISTMOTS'
      CALL LIROBJ(TYPE,ILIINC,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      MLMOTS = ILIINC
      SEGACT MLMOTS
      NC = MLMOTS.MOTS(/2)
      SEGDES MLMOTS
      IF(NC .NE. (IDIM+2+NESP))THEN
         MOTERR(1:40) = 'LISTINCO = ???'
         WRITE(IOIMP,*) MOTERR
C
C******* Message d'erreur standard
C        21 2
C        Données incompatibles
C
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C
C**** Boundary condition
C
      IRET=0
      TYPE='MAILLAGE'
      CALL LIROBJ(TYPE,IFLIM,0,IRET)
      IF(IERR.NE.0)GOTO 9999
      IF(IRET .EQ. 0)THEN
         MELLIM = 0
      ELSE
         MELEME=IFLIM
         SEGACT MELEME
         ICACCA=MELEME.NUM(/2)
         IF(ICACCA .EQ. 0)THEN
            MELLIM = 0
         ELSE
            MELLIM = IFLIM
         ENDIF
         SEGDES MELEME
      ENDIF
C
C**** Bas Mach (AUSMPLM, RUSANOLM, ROELM, HLLCLM, AUSMPUP)
C
      IF((INDMET .EQ. 9)  .OR. (INDMET .EQ. 11) .OR.
     &   (INDMET .EQ. 13) .OR. (INDMET .EQ. 15) .OR.
     &   (INDMET .EQ. 16)) THEN
         TYPE = 'CHPOINT '
C
C******* Reference speed
C
         CALL LIROBJ(TYPE,IUINF,1,IRET)
         IF(IERR .NE. 0) GOTO 9999
         INDIC = 1
         NBCOMP = 1
         MOT = 'SCAL'
         CALL QUEPOI(IUINF, MELEMC, INDIC, NBCOMP, MOT)
         IF(IERR .NE. 0) GOTO 9999
C
C******* Minimal cutoff
C
         TYPE = 'CHPOINT '
         CALL LIROBJ(TYPE,IUPRI,1,IRET)
         IF(IERR .NE. 0) GOTO 9999
         INDIC = 1
         NBCOMP = 1
         MOT = 'SCAL'
         CALL QUEPOI(IUPRI, MELEMC, INDIC, NBCOMP, MOT)
         IF(IERR .NE. 0) GOTO 9999
C
      ELSE
         IUINF=0
         IUPRI=0
      ENDIF
C
C**** Creation des flux aux interfaces
C
      TYPE = 'CHPOINT '
      CALL KRCHP1(TYPE, MELEMF, ICHFLU, MLMOTS)
C
C**** Calcul des flux et du pas du temps.
C
      IF(IDIM .EQ. 2)THEN
         CALL KONFL1(LOGME,INDMET,
     &        IROF,IVITF,IPF,IGAMF,IFRMAF,
     &        ICHPSU,ICHPDI,
     &        MELEMC,MELEMF,MELEFE,MELLIM,
     &        IUINF,IUPRI,
     &        ICHFLU,
     &        DT,
     &        LOGNC,LOGAN,MESERR)
      ELSEIF(IDIM .EQ. 3)THEN
         CALL KONFL3(LOGME,INDMET,
     &        IROF,IVITF,IPF,IGAMF,IFRMAF,
     &        ICHPSU,ICHPDI,
     &        MELEMC,MELEMF,MELEFE,MELLIM,
     &        IUINF,IUPRI,
     &        ICHFLU,
     &        DT,
     &        LOGNC,LOGAN,MESERR)
      ENDIF
C
      IF(LOGAN)THEN
C
C******* Anomalie detectée
C
C
C******* Message d'erreur standard
C        -301 0
C        %m1:40
C
         MOTERR(1:40) =  MESERR(1:40)
         WRITE(IOIMP,*) MOTERR(1:40)
C
C******* Message d'erreur standard
C        5 3
C        Erreur anormale.contactez votre support
C
         CALL ERREUR(5)
         GOTO 9999
      ENDIF
      IF(LOGNC)THEN
C
C******* Message d'erreur standard
C        -301 0
C        %m1:40
C
         MOTERR(1:40) =  MESERR(1:40)
         WRITE(IOIMP,*) MOTERR(1:40)
C
C******* Message d'erreur standard
C        460 2
C        Pas de convergence dans les itérations internes
C
         CALL ERREUR(460)
         GOTO 9999
      ENDIF
C
C**** Calcul de residu (IFLU = 2)
C
      IF(IFLU .EQ. 2)THEN
         TYPE = 'CHPOINT '
         CALL KRCHP1(TYPE, MELEMC, ICHRES, MLMOTS)
C
         CALL KONRE1(MELEMC,MELEMF,MELEFE,ICHPVO,
     &        ICHFLU, ICHRES,
     &        LOGAN,MESERR)
         IF(LOGAN)THEN
C
C******* Anomalie detectée
C
C
C******* Message d'erreur standard
C        -301 0
C        %m1:40
C
            MOTERR(1:40) =  MESERR(1:40)
            WRITE(IOIMP,*) MOTERR(1:40)
C
C******* Message d'erreur standard
C        5 3
C        Erreur anormale.contactez votre support
C
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
      ELSE
         ICHRES = 0
      ENDIF
C
C**** Ecriture des resultats
C
      CALL ECRREE(DT)
      TYPE = 'CHPOINT '
      IF(ICHRES .NE. 0) CALL ECROBJ(TYPE,ICHRES)
      IF(ICHFLU .NE. 0) CALL ECROBJ(TYPE,ICHFLU)
 9999 CONTINUE
      RETURN
      END