Télécharger la plaquette Cast3M


fr en

Télécharger cli272.eso

Retour à la liste

Numérotation des lignes : 

cli272

  1. C CLI272    SOURCE    OF166741  24/12/13    21:15:41     12097          
  2.       SUBROUTINE CLI272(NSP,MELEMF,MELEMC,MELECB,MELEFC,MELRES,INORM,
  3.      &     ICHPVO,ICHPSU,LRECP,LRECV,
  4.      &     IROC,IVITC,IPC,IYC,ICHLIM,ILIINC,ILIINP,IJAC,IJACO)
  5. C************************************************************************
  6. C
  7. C PROJET            :  CASTEM 2000
  8. C
  9. C NOM               :  CLI112
  10. C
  11. C DESCRIPTION       :  Subroutine appellée par CLIM22
  12. C                      OPTION: 'INJE' -- Jacobian
  13. C
  14. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  15. C
  16. C AUTEUR            :  S. Kudriakov, DEN/DM2S/SFME/LTMF
  17. C
  18. C************************************************************************
  19. C
  20. C APPELES (Calcul)  :
  21. C
  22. C************************************************************************
  23. C
  24. C HISTORIQUE (Anomalies et modifications éventuelles)
  25. C
  26. C HISTORIQUE :
  27. C
  28. C************************************************************************
  29. C
  30.       IMPLICIT INTEGER(I-N)
  31. C
  32.       INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
  33.      &     ,IGAMC,ICHLIM,ICEL,NFAC,IFAC,MELRES,IJACO
  34.      &     ,NGF,NGC,NLF,NLC,NLCB
  35.      &     ,ILIINC,ILIINP,IJAC,II,JJ,K
  36.      &     ,MP, NBEL, NBME, NBSOUS, NKID, NKMT, NMATRI, NP, NRIGE
  37.      &     ,NSP,I, IYC,J, LRECP,LRECV,KV
  38.       REAL*8 VOLU,SURF,RC,PC,UXC,UYC,GAMC,CNX,CNY
  39.      &     ,PSRF,RHOUF,DECIDP,DUNDPC,DUXDPC,DUYDPC
  40.      &     ,UN,RHO,P,COEF
  41.      &     ,BR1,BOT,TOP
  42.  
  43.       CHARACTER*(8) TYPE
  44.  
  45. -INC PPARAM
  46. -INC CCOPTIO
  47. -INC SMLMOTS
  48. -INC SMELEME
  49.       POINTEUR MELEFC.MELEME
  50. -INC SMLENTI
  51.       POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI
  52. -INC SMCHPOI
  53.       POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL,
  54.      &     MPVC.MPOVAL, MPPC.MPOVAL, MPLIM.MPOVAL, MPYC.MPOVAL
  56.       POINTEUR CELL.IZAFM
  57. C-------------------------------------------------------
  58. -INC SMLREEL
  59.       POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  60. C-------------------------------------------------------
  61. C*********  Les Jacobians ******************************
  62. C-------------------------------------------------------
  63.       SEGMENT JACEL
  64.          REAL*8 JAC(3+NSP,3+NSP)
  65.       ENDSEGMENT
  66.       POINTEUR JTL.JACEL, WL.JACEL, JTT.JACEL
  67. C-------------------------------------------------------------
  68. C******* Les fractionines massiques **************************
  69. C-------------------------------------------------------------
  70.        SEGMENT FRAMAS
  71.          REAL*8 YET(NSP)
  72.       ENDSEGMENT
  73.       POINTEUR YC.FRAMAS, YF.FRAMAS
  74. C-------------------------------------------------------
  75. C**********  Les CP's and CV's   ***********************
  76. C-------------------------------------------------------
  77.       SEGMENT GCONST
  78.          REAL*8 GC(NSP)
  79.       ENDSEGMENT
  80.       POINTEUR CP.GCONST, CV.GCONST
  81. C-------------------------------------------------------------
  82. C**********  Segments for the vectors  ***********************
  83. C-------------------------------------------------------------
  84.        SEGMENT VECEL
  85.          REAL*8 VV(NSP)
  86.        ENDSEGMENT
  87.        POINTEUR DGDYC.VECEL
  88. C--------------------------------------------------
  89.       SEGINI JTL
  90.       SEGINI JTT
  91.       SEGINI WL
  92. C----------------------------------------------------
  93. C**** KRIPAD pour la correspondance global/local
  94. C----------------------------------------------------
  95.       CALL KRIPAD(MELEMC,MLEMC)
  96.       CALL KRIPAD(MELECB,MLEMCB)
  97.       CALL KRIPAD(MELEMF,MLEMF)
  98. C----------------------------------------------------
  99. C**** CHPOINTs de la table DOMAINE
  100. C----------------------------------------------------
  101.       CALL LICHT(INORM,MPNORM,TYPE,ICEL)
  102.       CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL)
  103.       CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL)
  104. C----------------------------------------------------
  105. C**** CHPOINTs des variables
  106. C----------------------------------------------------
  107.       CALL LICHT(IROC,MPRC,TYPE,ICEL)
  108.       CALL LICHT(IVITC,MPVC,TYPE,ICEL)
  109.       CALL LICHT(IPC,MPPC,TYPE,ICEL)
  110.       CALL LICHT(IYC,MPYC,TYPE,ICEL)
  111.       CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
  112. C--------------------------------------------------------
  113. C**** Boucle sur le face pour le calcul des invariants de
  114. C     Riemann et du flux
  115. C--------------------------------------------------------
  116.       SEGACT MELEFC
  117.       NFAC=MELEFC.NUM(/2)
  118. C---------------------------------
  119. C**** Objet MATRIK
  120. C---------------------------------
  121.       NRIGE = 7
  122.       NMATRI = 1
  123.       NKID = 9
  124.       NKMT = 7
  125. C---------------------------------
  126.       SEGINI MATRIK
  127.       IJACO = MATRIK
  128.       MATRIK.IRIGEL(1,1) = MELRES
  129.       MATRIK.IRIGEL(2,1) = MELRES
  130. C---------------------------------
  131. C**** Matrice non symetrique
  132. C---------------------------------
  133.       MATRIK.IRIGEL(7,1) = 2
  134. C---------------------------------
  135.       NBME = (3+NSP)*(3+NSP)
  136.       NBSOUS = 1
  137.       SEGINI IMATRI
  138.       IF(IJAC.EQ.1)THEN
  139.          MLMOTS=ILIINC
  140.       ELSEIF(IJAC.EQ.2)THEN
  141.          MLMOTS=ILIINP
  142.       ENDIF
  143.       SEGACT MLMOTS
  144.       MATRIK.IRIGEL(4,1) = IMATRI
  145. C-------------------------------------------
  146.       DO 1 J=1,(NSP+3)
  147.         KV=(J-1)*(3+NSP)
  148.         IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1)
  149.         IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2)
  150.         IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3)
  151.         IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4)
  152.         DO 2 I=1,(NSP-1)
  153.           IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I)
  154.  2      CONTINUE
  155.  1    CONTINUE
  156. C-----------------------------------------------
  157.       SEGDES MLMOTS
  158.       MLMOTS=ILIINC
  159.       SEGACT MLMOTS
  160. C-----------------------------------------------
  161.       DO 3 J=1,(NSP+3)
  162.         KV=(J-1)*(3+NSP)
  163.         IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j)
  164.         IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j)
  165.         IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j)
  166.         IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j)
  167.         DO 4 I=1,(NSP-1)
  168.           IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j)
  169.  4     CONTINUE
  170.  3    CONTINUE
  171. C-----------------------------------------------
  172. C-----------------------------------------------
  173.       SEGDES MLMOTS
  174.       NBEL = NFAC
  175.       NBSOUS = 1
  176.       NP = 1
  177.       MP = 1
  178. C-----------------------------------------------------------
  179. C-----------------------------------------------------------
  180.       DO 5 I=1,NBME
  181.         SEGINI CELL
  182.         IMATRI.LIZAFM(1,I) = CELL
  183.  5    CONTINUE
  184. C---------------------------------
  185. C---------------------------------
  186. C**** Fin definition MATRIK
  187. C---------------------------------
  188.       DO IFAC=1,NFAC,1
  189.          NGF=MELEFC.NUM(1,IFAC)
  190.          NGC=MELEFC.NUM(2,IFAC)
  191.          NLF=MLEMF.LECT(NGF)
  192.          NLC=MLEMC.LECT(NGC)
  193.          NLCB=MLEMCB.LECT(NGF)
  194.          VOLU=MPVOL.VPOCHA(NLC,1)
  195.          SURF=MPSURF.VPOCHA(NLF,1)
  196. C        In CASTEM les normales sont sortantes
  197.          CNX=-1*MPNORM.VPOCHA(NLF,1)
  198.          CNY=-1*MPNORM.VPOCHA(NLF,2)
  199. C----------------------------------------------
  200.          SEGINI CP, CV
  201.          MLRECP = LRECP
  202.          MLRECV = LRECV
  203.          SEGACT MLRECP, MLRECV
  204.          DO 10 I=1,(NSP-1)
  205.            CP.GC(I)=MLRECP.PROG(I)
  206.            CV.GC(I)=MLRECV.PROG(I)
  207.  10      CONTINUE
  208.          CP.GC(NSP)=MLRECP.PROG(NSP)
  209.          CV.GC(NSP)=MLRECV.PROG(NSP)
  210. C---------------------------------
  211. C        Variables au centre
  212. C---------------------------------
  213.          RC=MPRC.VPOCHA(NLC,1)
  214.          PC=MPPC.VPOCHA(NLC,1)
  215.          UXC=MPVC.VPOCHA(NLC,1)
  216.          UYC=MPVC.VPOCHA(NLC,2)
  217.          SEGINI YC
  218.          SEGACT MPYC
  219.          DO 100 I=1,(NSP-1)
  220.            YC.YET(I)=MPYC.VPOCHA(NLC,I)
  221.  100     CONTINUE
  222. C---------------------------------
  223. C        Variables à la face
  224. C---------------------------------
  225.          RHOUF=MPLIM.VPOCHA(NLCB,1)
  226.          PSRF=MPLIM.VPOCHA(NLCB,2)
  227.          SEGINI YF
  228.          SEGACT MPLIM
  229.          DO 101 I=1,(NSP-1)
  230.            YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I)
  231.  101     CONTINUE
  232. c-------------------------------------------------------------
  233. c  Computing GAMMA at the cell-center
  234. c-------------------------------------------------------------
  235.          top=0.0D0
  236.          bot=0.0D0
  237.          do 102  i=1,(nsp-1)
  238.           top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
  239.           bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
  240.  102     continue
  241.          top=cp.gc(nsp)+top
  242.          bot=cv.gc(nsp)+bot
  243.          GAMC=top/bot
  244. C-------------------------------------------------------------
  245.          SEGINI DGDYC
  246.          do 41 i=1,(nsp-1)
  247.            dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  248.      &       GAMC*(cv.gc(i)-cv.gc(nsp)))/bot
  249.  41    continue
  250. c-------------------------------------------------------------
  251. C------------------------------------------------
  252. C******* Densite, vitesse, pression sur le bord
  253. C------------------------------------------------
  254.          P=PC
  255.          RHO=P/PSRF
  256.          UN=RHOUF/RHO
  257. C------------------------------
  258. C******* Derivatives w.r.t. PC
  259. C------------------------------
  260.          DUNDPC=-1*RHOUF/(RHO*RHO*PSRF)
  261.          DUXDPC=CNX*DUNDPC
  262.          DUYDPC=CNY*DUNDPC
  263.          DECIDP=UN*DUNDPC
  264. C--------------------------------------------------------------
  265.          COEF=SURF/VOLU
  266. C-------------------------------------------------------------
  267.          JTL.JAC(1,1) = 0.0D0
  268.          JTL.JAC(1,2) = 0.0D0
  269.          JTL.JAC(1,3) = 0.0D0
  270.          JTL.JAC(1,4) = 0.0D0
  271.          DO 107 I=1,(NSP-1)
  272.            JTL.JAC(1,4+I) = 0.0D0
  273.  107     CONTINUE
  274. C----------------------------------------
  275.          JTL.JAC(2,1) = 0.0D0
  276.          JTL.JAC(2,2) = 0.0D0
  277.          JTL.JAC(2,3) = 0.0D0
  278.          JTL.JAC(2,4) = (RHOUF*DUXDPC+CNX)*COEF
  279.          DO 108 I=1,(NSP-1)
  280.            JTL.JAC(2,4+I) = 0.0D0
  281.  108     CONTINUE
  282. C----------------------------------------
  283.          JTL.JAC(3,1) = 0.0D0
  284.          JTL.JAC(3,2) = 0.0D0
  285.          JTL.JAC(3,3) = 0.0D0
  286.          JTL.JAC(3,4) = (RHOUF*DUYDPC+CNY)*COEF
  287.          DO 109 I=1,(NSP-1)
  288.            JTL.JAC(3,4+I) = 0.0D0
  289.  109     CONTINUE
  290. C----------------------------------------
  291.          JTL.JAC(4,1) = 0.0D0
  292.          JTL.JAC(4,2) = 0.0D0
  293.          JTL.JAC(4,3) = 0.0D0
  294.          JTL.JAC(4,4) = (RHOUF*DECIDP)*COEF
  295.          DO 110 I=1,(NSP-1)
  296.            JTL.JAC(4,4+I) = 0.0D0
  297.  110     CONTINUE
  298. C----------------------------------------
  299.          DO 111 I=1,(NSP-1)
  300.           JTL.JAC(4+I,1) = 0.0D0
  301.           JTL.JAC(4+I,2) = 0.0D0
  302.           JTL.JAC(4+I,3) = 0.0D0
  303.           JTL.JAC(4+I,4) = 0.0D0
  304.           DO 112 J=1,(NSP-1)
  305.            JTL.JAC(4+I,4+J) = 0.0D0
  306.  112      CONTINUE
  307.  111     CONTINUE
  308. C---------------------------------------------------------
  309. c  matrix wl(i,j) represents the derivative of the i-component
  310. c  of the vector of primitive variables of the left state with
  311. c  respect to the j-component of the vector of the conservative
  312. c  variables of the left state.
  313. c
  314. c  Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  315. c                               vector of primitive variables;
  316. c  (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  317. c                            vector of conservative variables.
  318. c-------------------------------------------------------------
  319.        wl.jac(1,1)=1.0d0
  320.        wl.jac(1,2)=0.0d0
  321.        wl.jac(1,3)=0.0d0
  322.        wl.jac(1,4)=0.0d0
  323.        do 83 i=1,(nsp-1)
  324.          wl.jac(1,4+i)=0.0d0
  325.  83    continue
  326. c------------------------------
  327.        wl.jac(2,1)=-UXC/RC
  328.        wl.jac(2,2)=1.0d0/RC
  329.        wl.jac(2,3)=0.0d0
  330.        wl.jac(2,4)=0.0d0
  331.        do 84 i=1,(nsp-1)
  332.          wl.jac(2,4+i)=0.0d0
  333.  84    continue
  334. c------------------------------
  335.        wl.jac(3,1)=-UYC/RC
  336.        wl.jac(3,2)=0.0d0
  337.        wl.jac(3,3)=1.0d0/RC
  338.        wl.jac(3,4)=0.0d0
  339.        do 85 i=1,(nsp-1)
  340.          wl.jac(3,4+i)=0.0d0
  341.  85    continue
  342. c------------------------------
  343.        br1=0.0d0
  344.        do 86 i=1,(nsp-1)
  345.          br1=br1+dgdyc.vv(i)*yc.yet(i)
  346.  86    continue
  347.        br1=br1*PC/(RC*(GAMC-1.0D0))
  348.        wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1
  349.        wl.jac(4,2)=-UXC*(GAMC-1.0D0)
  350.        wl.jac(4,3)=-UYC*(GAMC-1.0D0)
  351.        wl.jac(4,4)=(GAMC-1.0D0)
  352.        do 87 i=1,(nsp-1)
  353.          wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0))
  354.  87    continue
  355. c------------------------------
  356.        do 88 i=1,(nsp-1)
  357.         do 89 j=1,4
  358.           wl.jac(4+i,j)=0.0d0
  359.           if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC
  360.  89     continue
  361. c------------
  362.         do 890 j=5,(4+nsp-1)
  363.           wl.jac(4+i,j)=0.0d0
  364.           if(4+i.eq.j)  then
  365.             wl.jac(4+i,j)=1.0d0/RC
  366.           endif
  367.  890    continue
  368.  88    continue
  369. c------------------------------------------------
  370. C------------------------------------------------
  371.        do 114 i=1,(3+nsp)
  372.         do 115 j=1,(3+nsp)
  373.          jtt.jac(i,j)=0.0d0
  374.          do 116 k=1,(3+nsp)
  375.           jtt.jac(i,j)=jtt.jac(i,j)+jtl.jac(i,k)*wl.jac(k,j)
  376.  116     continue
  377.  115    continue
  378.  114  continue
  379. C----------------------------------------------------------------
  380. C******* Jacobian with respect to conservative variables
  381. C----------------------------------------------------------------
  382.          IF(IJAC.EQ.1)THEN
  383.            DO 9 II = 1,(3+NSP)
  384.               DO 15 JJ = 1,(3+NSP)
  385.                 KV = (II-1)*(3+NSP)
  386. C----------------------------------
  387.                 CELL = IMATRI.LIZAFM(1,KV+JJ)
  388.                 CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ)
  389.  15           CONTINUE
  390.  9         CONTINUE
  391.          ELSEIF(IJAC.EQ.2)THEN
  392.            DO 20 II = 1,(3+NSP)
  393.               DO 25 JJ = 1,(3+NSP)
  394.                 KV = (II-1)*(3+NSP)
  395. C----------------------------------
  396.                 CELL = IMATRI.LIZAFM(1,KV+JJ)
  397.                 CELL.AM(IFAC,1,1) = JTL.JAC(II,JJ)
  398.  25           CONTINUE
  399.  20         CONTINUE
  400.          ENDIF
  401. c--------------------------------------------------
  402.       ENDDO
  403. C
  404.       SEGDES MELEFC
  405. C
  406.       SEGSUP MLEMC
  407.       SEGSUP MLEMCB
  408.       SEGSUP MLEMF
  409. C
  410.       SEGDES MPNORM
  411.       SEGDES MPVOL
  412.       SEGDES MPSURF
  413.       SEGDES MPRC
  414.       SEGDES MPPC
  415.       SEGDES MPVC
  416.       SEGDES MPYC
  417.       SEGDES MPLIM
  418.       SEGDES YC
  419.       SEGDES YF
  420.       SEGDES CP
  421.       SEGDES CV
  422.       SEGDES JTL
  423.       SEGDES JTT
  424.       SEGDES WL
  425.       SEGDES  DGDYC
  426.       SEGDES MATRIK
  427.       DO 80 II=1,NBME
  428.        CELL = IMATRI.LIZAFM(1,II)
  429.        SEGDES CELL
  430.  80   CONTINUE
  431.       SEGDES IMATRI
  432. C---------------------------------------------
  433.  9999 CONTINUE
  434.       RETURN
  435.       END
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  

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