Sources Esope | Cast3M

kbmms2
C KBMMS2    SOURCE    CB215821  20/11/25    13:30:48     10792          
      SUBROUTINE KBMMS2(NSP,ILIINC,IRN,IUN,IPN,IYN,MLRECP,
     &    MLRECV,INORM,ICHPVO,ICHPSU,IUINF,IUPRI,
     &    MELEMC,MELEFE,MELLIM,IMAT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KONMSP ("convection for multispecies")
C
C DESCRIPTION       :  Voir KON19 (appele par KON19)
C                      Calcul du jacobien du résidu pour la méthode
C                      AUSMplus
C
C                      Cas deux dimensions, gaz "calorically perfect"
C                                           MULTISPECIES!!!!!!!!
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  S. KUDRIAKOV, DM2S/SFME/LTMF
C
C************************************************************************
C
C
C APPELES (Outils
C          CASTEM)  :  KRIPAD, LICHT, ERREUR
C
C APPELES (Calcul)  :  CONMSP, CWMSP
C
C************************************************************************
C
C ENTREES
C
C     NSP    : number of species (total) ;
C
C     ILINC  : liste des inconnues (pointeur d'un objet de type LISTMOTS)
C
C  1) Pointeurs des CHPOINT
C
C     IRN     : CHPOINT CENTRE  contenant la masse volumique ;
C
C     IUN     : CHPOINT CENTRE  contenant la vitesse ;
C
C     IPN     : CHPOINT CENTRE  contenant la pression ;
C
C     IYN     : CHPOINT CENTRE  contenant les fractions massiques ;
C
C     INORM   : CHPOINT FACE    contenant les normales aux faces ;
C
C     ICHPOVO : CHPOINT VOLUME  contenant le volume
C
C     ICHPOSU : CHPOINT FACE    contenant la surface des faces
C
C
C  2) Pointeurs des LIST REELS
C
C     MLRECP  : LIST REELS contenant les CP's des gases differents ;
C
C     MLRECV  : LIST REELS contenant les CV's des gases differents ;
C
C  3) Pointeurs de MELEME de la table DOMAINE
C
C     MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
C
C     MELEFE  : MELEME 'FACEL' du connectivité Faces -> Elts
C
C     MELLIM  : MELEME SPG des conditions aux bords
C               N.B.: Sur le bord on ne calcule rien!
C
C SORTIES
C
C     IMAT    : pointeur de la MATRIK du jacobien du residu
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         GAMMA \in (1,3)
C       Si non il faut le faire!!!
C
C************************************************************************
C
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
C---------------------------------------------------
      INTEGER NSP,II,JJ,KV,I,J,JLL,JRR
      INTEGER ILIINC, IRN,IUN,IPN,IYN,INORM,ICHPVO,ICHPSU
     &        , IMAT, IGEOMC, IGEOMF, IUPRI, IUINF
     &        , NFAC, NBSOUS, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID
     &        , NKMT, NBME, NBEL,  MP, NP
     &        , IFAC, NGCF, NLCF, NGCG, NGCD, NLCG, NLCD, NLFL
      REAL*8 ROG, PG, UXG, UYG, VOLG
     &       , ROD, PD, UXD, UYD, VOLD
     &       , SURF, FUNCEL, UPR_L, UPR_R, UPR_M, V_INF
      REAL*8 WVEC_L(4), WVEC_R(4), NVECT(2), TVECT(2)
      REAL*8 ZC1, ZC2, ZC3, ZC4
      CHARACTER*8 TYPE
C
C**** LES INCLUDES
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMMATRIK
-INC SMELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MPRN.MPOVAL, MPUN.MPOVAL, MPPN.MPOVAL, MPYN.MPOVAL,
     &         MPNORM.MPOVAL, MPVOLU.MPOVAL, MPOVSU.MPOVAL
      POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME,
     &         MELEDU.MELEME, MELLIM.MELEME
      POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLELIM.MLENTI
      POINTEUR CELL.IZAFM
      POINTEUR MLMINC.MLMOTS, MPUPRI.MPOVAL, MPUINF.MPOVAL
C----------------------------------------------------
-INC SMLREEL
      POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
C-------------------------------------------------------
C*********  Les Jacobians ******************************
C-------------------------------------------------------
      SEGMENT JACEL
         REAL*8 JAC(3+NSP,3+NSP)
      ENDSEGMENT
      POINTEUR JTL.JACEL, JTR.JACEL, JTT.JACEL
C--------------------------------------------------------
C KRIPAD pour la correspondance global/local des centres
C--------------------------------------------------------
      CALL KRIPAD(MELLIM,MLELIM)
      CALL KRIPAD(MELEMC,MLENTC)
C------------------
      SEGACT MELEMC
C------------------
      SEGACT MELEFE
C------------------
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOMF)
      CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC)
C----------------------------------------------
      MELEMF = IGEOMF
      CALL KRIPAD(MELEMF,MLENTF)
      SEGACT MELEMF
C-----------------------------------------------
      CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
      CALL LICHT(IPN,MPPN,TYPE,IGEOMC)
      CALL LICHT(IUN,MPUN,TYPE,IGEOMC)
      CALL LICHT(IYN,MPYN,TYPE,IGEOMC)
C-----------------------------------------------
      NFAC = MELEFE.NUM(/2)
C-----------------------------------------------
C**** Maillage des inconnues primales
C-----------------------------------------------
      NBSOUS = 0
      NBREF = 0
      NBELEM =  NFAC
      NBNN = 2
C-----------------------------------------------
      SEGINI MELEDU
C----------------------
C**  MELEDU = 'SEG2' **
C----------------------
      MELEDU.ITYPEL =  2
C--------------
      NRIGE = 7
      NMATRI = 1
      NKID = 9
      NKMT = 7
C--------------
      SEGINI MATRIK
      IMAT = MATRIK
      MATRIK.IRIGEL(1,1) = MELEDU
      MATRIK.IRIGEL(2,1) = MELEDU
C-------------------------
C Matrice non symetrique
C-------------------------
      MATRIK.IRIGEL(7,1) = 2
C-------------------------
      NBME = (3+NSP)*(3+NSP)
      NBSOUS = 1
      SEGINI IMATRI
      MLMINC = ILIINC
      SEGACT MLMINC
      MATRIK.IRIGEL(4,1) = IMATRI
C-----------------------------------------------
      DO 1 J=1,(NSP+3)
        KV=(J-1)*(3+NSP)
        IMATRI.LISPRI(KV+1) = MLMINC.MOTS(1)
        IMATRI.LISPRI(KV+2) = MLMINC.MOTS(2)
        IMATRI.LISPRI(KV+3) = MLMINC.MOTS(3)
        IMATRI.LISPRI(KV+4) = MLMINC.MOTS(4)
        DO 2 I=1,(NSP-1)
          IMATRI.LISPRI(KV+4+I) = MLMINC.MOTS(4+I)
 2      CONTINUE
 1    CONTINUE
C-----------------------------------------------
      DO 3 J=1,(NSP+3)
        KV=(J-1)*(3+NSP)
        IMATRI.LISDUA(KV+1) = MLMINC.MOTS(j)
        IMATRI.LISDUA(KV+2) = MLMINC.MOTS(j)
        IMATRI.LISDUA(KV+3) = MLMINC.MOTS(j)
        IMATRI.LISDUA(KV+4) = MLMINC.MOTS(j)
        DO 4 I=1,(NSP-1)
          IMATRI.LISDUA(KV+4+I) = MLMINC.MOTS(j)
 4     CONTINUE
 3    CONTINUE
C-----------------------------------------------
C-----------------------------------------------
      NBEL = NBELEM
      NBSOUS = 1
      NP = 2
      MP = 2
C-----------------------------------------------------------
C-----------------------------------------------------------
      DO 5 I=1,NBME
        SEGINI CELL
        IMATRI.LIZAFM(1,I) = CELL
 5    CONTINUE
C**************************************************************
C    Bas Mach
C**************************************************************
      CALL LICHT(IUPRI,MPUPRI,TYPE,IGEOMC)
      CALL LICHT(IUINF,MPUINF,TYPE,IGEOMC)
C---------------------------------
      DO IFAC = 1, NFAC, 1
         NGCF = MELEFE.NUM(2,IFAC)
         NLCF = MLENTF.LECT(NGCF)
         IF(NLCF .NE. IFAC)THEN
            WRITE(IOIMP,*) 'Il ne faut pas jouer avec la table domaine'
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
         NLFL = MLELIM.LECT(NGCF)
         NGCG = MELEFE.NUM(1,IFAC)
         NGCD = MELEFE.NUM(3,IFAC)
         IF(NLFL .NE. 0)THEN
C---------------------------------------------------------
C The point belongs to BC -> No contribution to jacobian!
C---------------------------------------------------------
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
         ELSEIF(NGCG .NE. NGCD)THEN
C-----------------------------------
C********** Les MELEMEs
C-----------------------------------
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
C-----------------------------------
C********** Les etats G et D
C-----------------------------------
            NLCG = MLENTC.LECT(NGCG)
            NLCD = MLENTC.LECT(NGCD)
C-----------------
            ROG = MPRN.VPOCHA(NLCG,1)
            PG = MPPN.VPOCHA(NLCG,1)
            UXG = MPUN.VPOCHA(NLCG,1)
            UYG = MPUN.VPOCHA(NLCG,2)
            VOLG = MPVOLU.VPOCHA(NLCG,1)
C-------------------------------------------------
            WVEC_L(1)=ROG
            WVEC_L(2)=UXG
            WVEC_L(3)=UYG
            WVEC_L(4)=PG
C-------------------------------------------------
            ROD = MPRN.VPOCHA(NLCD,1)
            PD = MPPN.VPOCHA(NLCD,1)
            UXD = MPUN.VPOCHA(NLCD,1)
            UYD = MPUN.VPOCHA(NLCD,2)
            VOLD = MPVOLU.VPOCHA(NLCD,1)
C------------------------------------------------
            UPR_L=MPUPRI.VPOCHA(NLCG,1)
            UPR_R=MPUPRI.VPOCHA(NLCD,1)
            UPR_M=MAX(UPR_L,UPR_R)
            V_INF=MAX(MPUINF.VPOCHA(NLCG,1),MPUINF.VPOCHA(NLCD,1))
            V_INF=MAX(UPR_M,V_INF)
C------------------------------------------------
            WVEC_R(1)=ROD
            WVEC_R(2)=UXD
            WVEC_R(3)=UYD
            WVEC_R(4)=PD
C------------------------------------------------
C  La normale G->D
C  La tangente
C------------------------------------------------
            SURF = MPOVSU.VPOCHA(NLCF,1)
            NVECT(1) = MPNORM.VPOCHA(NLCF,1)
            NVECT(2) = MPNORM.VPOCHA(NLCF,2)
            TVECT(1) = -1.0D0 * NVECT(2)
            TVECT(2) = NVECT(1)
C-------------------------------------------------
            CALL jbmms2(NSP,JLL,JRR,WVEC_L,WVEC_R,NVECT,TVECT,
     &        MPYN,MLRECP,MLRECV,NLCG,NLCD,V_INF)
C-----------------------------------------------------
C********** AB.AM(IFAC,IPRIM,IDUAL)
C           A = nom de l'inconnu duale (Ro,rUX,rUY,RET)
C           B = nom de l'inconnu primale (Ro,rUX,rUY,RET)
C           IPRIM = 1, 2  -> G, D
C           IDUAL = 1, 2  -> G, D
C           i.e.
C           A_IDUAL = AB.AM(IFAC,IPRIM,IDUAL) * B_IPRIM + ...
C
C
C********** Dual RN
C----------------------------------------------------
            JTL=JLL
            JTR=JRR
            SEGACT JTL
            SEGACT JTR
            DO 10 II = 1,(3+NSP)
              DO 15 JJ = 1,(3+NSP)
                KV = (II-1)*(3+NSP)
C----------------------------------
                CELL = IMATRI.LIZAFM(1,KV+JJ)
                FUNCEL = SURF * JTL.JAC(II,JJ)
                CELL.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
                CELL.AM(IFAC,1,2) = FUNCEL / VOLD
 15           CONTINUE
 10         CONTINUE
C-----------------------------------------------------
            DO 20 II = 1,(3+NSP)
              DO 25 JJ = 1,(3+NSP)
                KV = (II-1)*(3+NSP)
                CELL = IMATRI.LIZAFM(1,KV+JJ)
                FUNCEL = SURF * JTR.JAC(II,JJ)
                CELL.AM(IFAC,2,2) = FUNCEL / VOLD
                CELL.AM(IFAC,2,1) = -FUNCEL / VOLG
 25           CONTINUE
 20         CONTINUE
            SEGDES JTR
            SEGDES JTL
C-----------------------------------------------------
         ELSE
C-----------------------------------------------------
C*************  Murs (NGCG = NGCD)  ******************
C-----------------------------------------------------
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
            NLCG = MLENTC.LECT(NGCG)
C--------------------------------------
            ROG = MPRN.VPOCHA(NLCG,1)
            PG = MPPN.VPOCHA(NLCG,1)
            UXG = MPUN.VPOCHA(NLCG,1)
            UYG = MPUN.VPOCHA(NLCG,2)
            VOLG = MPVOLU.VPOCHA(NLCG,1)
C-------------------------------------------
            WVEC_L(1)=ROG
            WVEC_L(2)=UXG
            WVEC_L(3)=UYG
            WVEC_L(4)=PG
C-------------------------------------------------
            SURF = MPOVSU.VPOCHA(NLCF,1)
            NVECT(1) = MPNORM.VPOCHA(NLCF,1)
            NVECT(2) = MPNORM.VPOCHA(NLCF,2)
            TVECT(1) =-NVECT(2)
            TVECT(2) = NVECT(1)
C------- COEFFICIENTS ----------------------------
            ZC1=NVECT(1)*TVECT(2)+TVECT(1)*NVECT(2)
            ZC2=NVECT(1)*TVECT(2)-TVECT(1)*NVECT(2)
            ZC3=2.0D0*NVECT(1)*TVECT(1)
            ZC4=2.0D0*NVECT(2)*TVECT(2)
C-------------------------------------------------
            ROD = ROG
            PD = PG
            UXD = -ZC1*UXG/ZC2-ZC4*UYG/ZC2
            UYD = ZC3*UXG/ZC2+ZC1*UYG/ZC2
            VOLD = VOLG
C------------------------------------------------
            UPR_L=MPUPRI.VPOCHA(NLCG,1)
            UPR_R=UPR_L
            V_INF=MAX(MPUINF.VPOCHA(NLCG,1),UPR_L)
C------------------------------------------------
            WVEC_R(1)=ROD
            WVEC_R(2)=UXD
            WVEC_R(3)=UYD
            WVEC_R(4)=PD
C-------------------------------------------
C********** La normale sortante
C-------------------------------------------
             CALL jbmmw2(NSP,JLL,WVEC_L,WVEC_R,NVECT,
     &        TVECT,MPYN,MLRECP,MLRECV,NLCG,V_INF)
C----------------------------------------------------
            JTL=JLL
            SEGACT JTL
            DO 70 II = 1,(3+NSP)
              DO 75 JJ = 1,(3+NSP)
                KV = (II-1)*(3+NSP)
C---------------------------------
                CELL = IMATRI.LIZAFM(1,KV+JJ)
                FUNCEL = SURF * JTL.JAC(II,JJ)
                CELL.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
                CELL.AM(IFAC,1,2) = 0.0D0
 75           CONTINUE
 70         CONTINUE
            SEGDES JTL
C--------------------------------------------
            DO 40 II = 1,(3+NSP)
              DO 45 JJ = 1,(3+NSP)
                KV = (II-1)*(3+NSP)
                CELL = IMATRI.LIZAFM(1,KV+JJ)
                CELL.AM(IFAC,2,2) = 0.0D0
                CELL.AM(IFAC,2,1) = 0.0D0
 45           CONTINUE
 40         CONTINUE
C--------------------------------------------
         ENDIF
      ENDDO
C
      SEGDES MELEMC
      SEGDES MELEFE
      SEGDES MELEMF
C
      SEGDES MPOVSU
      SEGDES MPVOLU
      SEGDES MPNORM
C
      SEGDES MPRN
      SEGDES MPPN
      SEGDES MPUN
      SEGDES MPYN
C
      SEGDES MELEDU
      SEGDES MATRIK
      DO 80 II=1,NBME
       CELL = IMATRI.LIZAFM(1,II)
       SEGDES CELL
 80   CONTINUE
      SEGDES IMATRI
C
      SEGSUP MLENTC
      SEGSUP MLENTF
      SEGSUP MLELIM
      SEGDES MLMINC
C
      SEGDES MPUPRI
      SEGDES MPUINF
      IF(MELLIM .NE.0) SEGDES MELLIM
 
 9999 CONTINUE
      RETURN
      END