Télécharger kon13.eso

Retour à la liste

kon13
C KON13     SOURCE    OF166741  24/12/13    21:16:13     12097          
      SUBROUTINE KON13
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KON13
C
C DESCRIPTION       :  Subroutine appellée par KON1
C
C                      Modelisation 2D/3D des equations d'Euler
C                      Calcul du jacobien (2eme ordre)
C                      Option non documentée dans la notice!!!
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)  :  KONJA5 (calcul du jacobien, gaz "calorically
C                              perfect", monoespece, 2D, 2eme ordre,
C                              AUSMplus, VLH)
C
C************************************************************************
C
C*** SYNTAXE
C
C  Discrétisation en VF "cell-centered" des équations d'Euler pour
C  un gaz parfait mono-constituent polytropique
C  Inconnues: densité, quantité de mouvement, énergie totale par
C  unité de volumes (variables conservatives)
C  Calcul du jacobien 2eme ordre en espace
C  Cas 2D
C
C  RMAT1 = 'KONV' 'VF' 'PERFMONO' 'JACOCON2' LMOT1
C          MOD1 MOT1
C          MCHAM1 MCHAM2 MCHAM3 MCHAM4
C          (MAILIM) CHPO1 CHPO2 CHPO3 CHPO5
C          CHPO6  CHPO7  CHPO8  CHPO9  CHPO10 CHPO11
C          MCHAM5 MCHAM6 MCHAM7 CHPO12 ;
C
C  ENTREES
C
C  MOT1   : objet de type MOT
C           'VLH'       : jacobien du residu pour la methode VLH
C           'AUSMPLUS'  : jacobien du residu pour la methode AUSM+
C
C  LMOT1  : objet de type LISTMOTS
C           Noms de composantes des variable primales et duales de RMAT1.
C           Il contient dans l'ordre suivant: le noms de la densité,
C           du momentum, de l'énergie totale par unité de volume
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 (MAILIM): MAIILAGE de POI1 ou on ne veut pas calculer le FLUX convective
C
C  CHPO1   : CHPOINT contenant la masse volumique
C           (SPG =(DOMA MOD1 'CENTRE'), une seule composante,
C            'SCAL').
C
C  CHPO2   : CHPOINT contenant la vitesse
C           (SPG =(DOMA MOD1 'CENTRE'), deux/trois composantes
C            'UX', 'UY', 'UZ')
C
C  CHPO3   : CHPOINT contenant la pression du gaz
C           (SPG =(DOMA MOD1 'CENTRE'), une seule composante,
C            'SCAL').
C
C  CHPO5   : CHPOINT contenant le "gamma" du gaz
C           (SPG =(DOMA MOD1 'CENTRE'), une seule composante,
C            'SCAL').
C
C  CHPO6   : CHPOINT contenant le gradient de la densité
C           CHPO1 (voir operatéur PENT)
C
C  CHPO7   : CHPOINT contenant le gradient de la vitesse
C           CHPO2 (voir operatéur PENT)
C
C  CHPO8   : CHPOINT contenant le gradient de la pression
C           CHPO3 (voir operatéur PENT)
C
C  CHPO9   : CHPOINT contenant le limiteur de la densité
C           CHPO1 (voir operatéur PENT)
C
C  CHPO10  : CHPOINT contenant le limiteur de la vitesse
C           CHPO2 (voir operatéur PENT)
C
C  CHPO11  : CHPOINT contenant le limiteur de la pression
C           CHPO3 (voir operatéur PENT)
C
C  MCHAM5  : CHAMELEM contenant le coefficient pour le calcul
C           du gradient de la densité
C           CHPO1 (voir operatéur PENT)
C
C  MCHAM6  : CHAMELEM contenant le coefficient pour le calcul
C           du gradient de la vitesse
C           CHPO2 (voir operatéur PENT)
C
C  MCHAM7  : CHAMELEM contenant le coefficient pour le calcul
C           du gradient de la pression
C           CHPO3 (voir operatéur PENT)
C
C  CHPO12  : CHPOINT contenant la condition limite sur la vitesse
C            aux bords
C
C SORTIES
C
C
C RMAT1   : objet de type MATRIK
C           (SPG =(DOMA MOD1 'CENTRE'))
C           (inconnues primales = inconnues duales = LMOT1)
C***********************************************************************
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
 
-INC PPARAM
-INC CCOPTIO
-INC SMLMOTS
-INC SMCHPOI
-INC SMELEME
      POINTEUR MLMVIT.MLMOTS
C
      INTEGER NBJAC, IRET, INDIC, NBCOMP, NESP, NC, JGN, JGM
     &     ,IIMPL, IIMPL2, IJACO
     &     ,IDOMA, MELEMC, MELEMF, MELEFE, MELTFA, ICHPSU, ICHPDI
     &     ,ICHPVO, INORM
     &     ,IROF, IVITF, IPF, IGAMF, ILIINC
     &     ,IFLIM, ICACCA, MELLIM
     &     ,IRN, IVN, IPN, IGAMN
     &     ,IGRN, IGVN, IGPN, IGLRN, IGLVN, IGLPN, ICRN, ICVN, ICPN
     &     ,IVLIM, INEFMD, ICOND
C
      PARAMETER (NBJAC=2)
      CHARACTER*8 TYPE, LJACO(NBJAC)
      CHARACTER*4 MOT
C
      DATA LJACO/'VLH     ','AUSMPLUS'/
C
      CALL LIRMOT(LJACO,NBJAC,IIMPL,1)
      IF(IERR .NE. 0)GOTO 9999
C
C**********************************
C**** Lecture de l'objet MODELE ***
C**********************************
C
      ICOND = 1
      CALL QUETYP(TYPE,ICOND,IRET)
 
      IF(IRET.EQ.0.AND.TYPE.NE.'MMODEL')THEN
        WRITE(6,*)' On attend un objet MMODEL'
        RETURN
      ENDIF
      CALL LIROBJ('MMODEL',MMODEL,ICOND,IRET)
      IF(IERR.NE.0)GOTO 9999
      CALL LEKMOD(MMODEL,IDOMA,INEFMD)
      IF(IERR.NE.0)GOTO 9999
C
C**** Centre, FACE, FACEL, ELTFA
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
      CALL LEKTAB(IDOMA,'ELTFA',MELTFA)
      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
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**** La list des inconnues
C
      NESP=0
      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******* La densité au centre
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IRN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C**** Control du CHPOINT: QUEPOI
C
C     INDIC = 1   -> on impose le pointeur du support geometrique (ICEN)
C                   N.B. Le CHPOINT peut changer de structure pour
C                        avoir SPG = ICEN!!!!
C     INDIC = 0   -> on ne fait que verifier le support geometrique
C                   (ICEN). Si le SPG sont differents INDIC = -4 en sortie
C
C     NBCOMP > 0 -> numero des composantes
C
C     MOT  = ' ' obligatoire s'on connais pas les noms des composantes
C
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(IRN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C******* La vitesse au centre
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IVN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      JGN = 4
      JGM = IDIM
      SEGINI MLMVIT
      MLMVIT.MOTS(1) = 'UX  '
      MLMVIT.MOTS(2) = 'UY  '
      IF(IDIM .EQ. 3) MLMVIT.MOTS(3) = 'UZ  '
      CALL QUEPO1(IVN, MELEMC, MLMVIT)
      SEGSUP MLMVIT
      IF(IERR .NE. 0) GOTO 9999
C
C******* La pression au centre
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IPN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(IPN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C******* Gamma au centre
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGAMN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'SCAL'
      CALL QUEPOI(IGAMN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C
C********** Gradient of density
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGRN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      JGN = 4
      JGM = IDIM
      SEGINI MLMVIT
      MLMVIT.MOTS(1) = 'P1DX'
      MLMVIT.MOTS(2) = 'P1DY'
      IF(IDIM .EQ. 3) MLMVIT.MOTS(3) = 'P1DZ'
      CALL QUEPO1(IGRN, MELEMC, MLMVIT)
      SEGSUP MLMVIT
      IF(IERR .NE. 0) GOTO 9999
C
C********** Gradient of speed
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGVN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      JGN = 4
      JGM = IDIM*IDIM
      SEGINI MLMVIT
      MLMVIT.MOTS(1) = 'P1DX'
      MLMVIT.MOTS(2) = 'P1DY'
      MLMVIT.MOTS(3) = 'P2DX'
      MLMVIT.MOTS(4) = 'P2DY'
      IF(IDIM .EQ. 3)THEN
         MLMVIT.MOTS(1) = 'P1DX'
         MLMVIT.MOTS(2) = 'P1DY'
         MLMVIT.MOTS(3) = 'P1DZ'
         MLMVIT.MOTS(4) = 'P2DX'
         MLMVIT.MOTS(5) = 'P2DY'
         MLMVIT.MOTS(6) = 'P2DZ'
         MLMVIT.MOTS(7) = 'P3DX'
         MLMVIT.MOTS(8) = 'P3DY'
         MLMVIT.MOTS(9) = 'P3DZ'
      ENDIF
      CALL QUEPO1(IGVN, MELEMC, MLMVIT)
      SEGSUP MLMVIT
      IF(IERR .NE. 0) GOTO 9999
C
C********** Gradient of pressure
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGPN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      JGN = 4
      JGM = IDIM
      SEGINI MLMVIT
      MLMVIT.MOTS(1) = 'P1DX'
      MLMVIT.MOTS(2) = 'P1DY'
      IF(IDIM .EQ. 3) MLMVIT.MOTS(3) = 'P1DZ'
      CALL QUEPO1(IGPN, MELEMC, MLMVIT)
      SEGSUP MLMVIT
      IF(IERR .NE. 0) GOTO 9999
C
C********** Limiter of density gradient
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGLRN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'P1  '
      CALL QUEPOI(IGLRN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C********** Limiter of velocity gradient
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGLVN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      JGN = 4
      JGM = IDIM
      SEGINI MLMVIT
      MLMVIT.MOTS(1) = 'P1  '
      MLMVIT.MOTS(2) = 'P2  '
      IF(IDIM .EQ. 3) MLMVIT.MOTS(3) = 'P3  '
      CALL QUEPO1(IGLVN, MELEMC, MLMVIT)
      SEGSUP MLMVIT
      IF(IERR .NE. 0) GOTO 9999
C
C********** Limiter of pressure gradient
C
      TYPE = 'CHPOINT '
      CALL LIROBJ(TYPE,IGLPN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
      INDIC = 1
      NBCOMP = 1
      MOT = 'P1  '
      CALL QUEPOI(IGLPN, MELEMC, INDIC, NBCOMP, MOT)
      IF(IERR .NE. 0) GOTO 9999
C
C********** On va lire les pointeurs des MCHAMLs pour
C           le calcul du gradient
C
C           Densité
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,ICRN,1,IRET)
      IF(IERR.NE.0) GOTO 9999
C
C********** Vitesse
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,ICVN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C********** Pression
C
      TYPE='MCHAML  '
      CALL  LIROBJ(TYPE,ICPN,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C
C********** On va lire le chpoin limite pour la vitesse
C
      TYPE='CHPOINT '
      CALL LIROBJ(TYPE,IVLIM,1,IRET)
      IF(IERR .NE. 0) GOTO 9999
C********** CHPOINT vide?
      MCHPOI = IVLIM
      SEGACT MCHPOI
      ICACCA = MCHPOI.IPCHP(/1)
      SEGDES MCHPOI
      IF(ICACCA .EQ. 0)THEN
         IVLIM=0
      ELSE
         JGN = 4
         JGM = IDIM
         SEGINI MLMVIT
         MLMVIT.MOTS(1) = 'UX  '
         MLMVIT.MOTS(2) = 'UY  '
         IF(IDIM .EQ. 3) MLMVIT.MOTS(3) = 'UZ  '
         ICACCA=0
         CALL QUEPO1(IVLIM, ICACCA, MLMVIT)
         IF(IERR .NE. 0) GOTO 9999
         SEGSUP MLMVIT
      ENDIF
C
C******* Calcul du jacobien
C
      IF(IDIM .EQ. 2)THEN
C        IIMPL = 4 VLH
C        IIMPL = 5 AUSM+
         IIMPL2 = IIMPL+3
         CALL KONJA5(IIMPL2,ILIINC,IROF,IPF,IVITF,IGAMF,
     &        IRN,IPN,IVN,IGAMN,
     &        ICHPVO,ICHPSU, INORM, MELEMC, MELEFE, MELEMF, MELTFA
     &        ,IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN,ICRN, ICVN,
     &        ICPN, IVLIM, MELLIM, IJACO)
         IF(IERR .NE. 0) GOTO 9999
      ELSE
C              Tentative d'utilisation d'une option non implémentée
         CALL ERREUR(251)
         GOTO 9999
      ENDIF
C
      TYPE='MATRIK '
      CALL ECROBJ(TYPE,IJACO)
 9999 CONTINUE
      RETURN
      END
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales