Télécharger la plaquette Cast3M


fr en

Télécharger cli262.eso

Retour à la liste

Numérotation des lignes : 

cli262

  1. C CLI262    SOURCE    OF166741  24/12/13    21:15:38     12097          
  2.       SUBROUTINE CLI262(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: 'INSU' -- 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.       INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
  32.      &     ,IGAMC,ICHLIM,ICEL,NFAC,IFAC,MELRES,IJACO
  33.      &     ,NGF,NGC,NLF,NLC,NLCB
  34.      &     ,ILIINC,ILIINP,IJAC,II,JJ,K
  35.      &     ,MP, NBEL, NBME, NBSOUS, NKID, NKMT, NMATRI, NP, NRIGE
  36.      &     ,NSP,I, IYC,J, LRECP,LRECV,KV
  37.       REAL*8 VOLU,SURF,RC,PC,UXC,UYC,GAMC,CNX,CNY,CTX,CTY
  38.      &     ,PSRF,COEF
  39.      &     ,BR1,BOT,TOP,DUNDUY,GAMF,UNC
  40.      &     ,COEF5,DPSRUN,DUXDUN,DRDUN,DUYDUN,DUNDUX,DPDUN
  41.      &     ,GM1,USGM1,RF,UNF,UTF,UXF,UYF,HTF,PF,SF,ECIN
  42.      &     ,TVECT(2),NVECT(2),WVEC_L(4),WVEC_R(4)
  43.       CHARACTER*(8) TYPE
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC SMLMOTS
  47. -INC SMELEME
  48.       POINTEUR MELEFC.MELEME
  49. -INC SMLENTI
  50.       POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI
  51. -INC SMCHPOI
  52.       POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL,
  53.      &     MPVC.MPOVAL, MPPC.MPOVAL, MPLIM.MPOVAL, MPYC.MPOVAL
  55.       POINTEUR CELL.IZAFM
  56. C-------------------------------------------------------
  57. -INC SMLREEL
  58.       POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  59. C-------------------------------------------------------
  60. C*********  Les Jacobians ******************************
  61. C-------------------------------------------------------
  62.       SEGMENT JACEL
  63.          REAL*8 JAC(3+NSP,3+NSP)
  64.       ENDSEGMENT
  65.       POINTEUR JLL.JACEL, WL.JACEL, JRR.JACEL, JTT.JACEL,
  66.      &         JFL.JACEL, JFR.JACEL, JR.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 JTT
  90.       SEGINI JLL
  91.       SEGINI JRR
  92.       SEGINI JFL
  93.       SEGINI JFR
  94.       SEGINI WL
  95. C----------------------------------------------------
  96. C**** KRIPAD pour la correspondance global/local
  97. C----------------------------------------------------
  98.       CALL KRIPAD(MELEMC,MLEMC)
  99.       CALL KRIPAD(MELECB,MLEMCB)
  100.       CALL KRIPAD(MELEMF,MLEMF)
  101. C----------------------------------------------------
  102. C**** CHPOINTs de la table DOMAINE
  103. C----------------------------------------------------
  104.       CALL LICHT(INORM,MPNORM,TYPE,ICEL)
  105.       CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL)
  106.       CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL)
  107. C----------------------------------------------------
  108. C**** CHPOINTs des variables
  109. C----------------------------------------------------
  110.       CALL LICHT(IROC,MPRC,TYPE,ICEL)
  111.       CALL LICHT(IVITC,MPVC,TYPE,ICEL)
  112.       CALL LICHT(IPC,MPPC,TYPE,ICEL)
  113.       CALL LICHT(IYC,MPYC,TYPE,ICEL)
  114.       CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
  115. C--------------------------------------------------------
  116. C**** Boucle sur le face pour le calcul des invariants de
  117. C     Riemann et du flux
  118. C--------------------------------------------------------
  119.       SEGACT MELEFC
  120.       NFAC=MELEFC.NUM(/2)
  121. C---------------------------------
  122. C**** Objet MATRIK
  123. C---------------------------------
  124.       NRIGE = 7
  125.       NMATRI = 1
  126.       NKID = 9
  127.       NKMT = 7
  128. C---------------------------------
  129.       SEGINI MATRIK
  130.       IJACO = MATRIK
  131.       MATRIK.IRIGEL(1,1) = MELRES
  132.       MATRIK.IRIGEL(2,1) = MELRES
  133. C---------------------------------
  134. C**** Matrice non symetrique
  135. C---------------------------------
  136.       MATRIK.IRIGEL(7,1) = 2
  137. C---------------------------------
  138.       NBME = (3+NSP)*(3+NSP)
  139.       NBSOUS = 1
  140.       SEGINI IMATRI
  141.       IF(IJAC.EQ.1)THEN
  142.          MLMOTS=ILIINC
  143.       ELSEIF(IJAC.EQ.2)THEN
  144.          MLMOTS=ILIINP
  145.       ENDIF
  146.       SEGACT MLMOTS
  147.       MATRIK.IRIGEL(4,1) = IMATRI
  148. C-------------------------------------------
  149.       DO 1 J=1,(NSP+3)
  150.         KV=(J-1)*(3+NSP)
  151.         IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1)
  152.         IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2)
  153.         IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3)
  154.         IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4)
  155.         DO 2 I=1,(NSP-1)
  156.           IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I)
  157.  2      CONTINUE
  158.  1    CONTINUE
  159. C-----------------------------------------------
  160.       SEGDES MLMOTS
  161.       MLMOTS=ILIINC
  162.       SEGACT MLMOTS
  163. C-----------------------------------------------
  164.       DO 3 J=1,(NSP+3)
  165.         KV=(J-1)*(3+NSP)
  166.         IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j)
  167.         IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j)
  168.         IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j)
  169.         IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j)
  170.         DO 4 I=1,(NSP-1)
  171.           IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j)
  172.  4     CONTINUE
  173.  3    CONTINUE
  174. C-----------------------------------------------
  175. C-----------------------------------------------
  176.       SEGDES MLMOTS
  177.       NBEL = NFAC
  178.       NBSOUS = 1
  179.       NP = 1
  180.       MP = 1
  181. C-----------------------------------------------------------
  182. C-----------------------------------------------------------
  183.       DO 5 I=1,NBME
  184.         SEGINI CELL
  185.         IMATRI.LIZAFM(1,I) = CELL
  186.  5    CONTINUE
  187. C---------------------------------
  188. C---------------------------------
  189. C**** Fin definition MATRIK
  190. C---------------------------------
  191.       DO IFAC=1,NFAC,1
  192.          NGF=MELEFC.NUM(1,IFAC)
  193.          NGC=MELEFC.NUM(2,IFAC)
  194.          NLF=MLEMF.LECT(NGF)
  195.          NLC=MLEMC.LECT(NGC)
  196.          NLCB=MLEMCB.LECT(NGF)
  197.          VOLU=MPVOL.VPOCHA(NLC,1)
  198.          SURF=MPSURF.VPOCHA(NLF,1)
  199. C        In CASTEM les normales sont sortantes
  200.          CNX=-1*MPNORM.VPOCHA(NLF,1)
  201.          CNY=-1*MPNORM.VPOCHA(NLF,2)
  202.          CTX=-1.0D0*CNY
  203.          CTY=CNX
  204. C----------------------------------------------
  205.          SEGINI CP, CV
  206.          MLRECP = LRECP
  207.          MLRECV = LRECV
  208.          SEGACT MLRECP, MLRECV
  209.          DO 10 I=1,(NSP-1)
  210.            CP.GC(I)=MLRECP.PROG(I)
  211.            CV.GC(I)=MLRECV.PROG(I)
  212.  10      CONTINUE
  213.          CP.GC(NSP)=MLRECP.PROG(NSP)
  214.          CV.GC(NSP)=MLRECV.PROG(NSP)
  215. C---------------------------------
  216. C        Variables au centre
  217. C---------------------------------
  218.          RC=MPRC.VPOCHA(NLC,1)
  219.          PC=MPPC.VPOCHA(NLC,1)
  220.          UXC=MPVC.VPOCHA(NLC,1)
  221.          UYC=MPVC.VPOCHA(NLC,2)
  222.          SEGINI YC
  223.          SEGACT MPYC
  224.          DO 100 I=1,(NSP-1)
  225.            YC.YET(I)=MPYC.VPOCHA(NLC,I)
  226.  100     CONTINUE
  227.          UNC=(UXC*CNX)+(UYC*CNY)
  228. C---------------------------------
  229. C        Variables à la face
  230. C---------------------------------
  231.          SEGACT MPLIM
  232.          HTF=MPLIM.VPOCHA(NLCB,1)
  233.          SF=MPLIM.VPOCHA(NLCB,2)
  234.          SEGINI YF
  235.          DO 101 I=1,(NSP-1)
  236.            YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I)
  237.  101     CONTINUE
  238.          UNF=UNC
  239.          UTF=0.0D0
  240.          UXF=UNF*CNX+UTF*CTX
  241.          UYF=UNF*CNY+UTF*CTY
  242. c-------------------------------------------------------------
  243. c  Computing GAMMA at the face
  244. c-------------------------------------------------------------
  245.          top=0.0D0
  246.          bot=0.0D0
  247.          do 102  i=1,(nsp-1)
  248.           top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp))
  249.           bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp))
  250.  102     continue
  251.          top=cp.gc(nsp)+top
  252.          bot=cv.gc(nsp)+bot
  253.          GAMF=top/bot
  254.          GM1=GAMF-1.0D0
  255.          USGM1=1.0D0/GM1
  256. c-------------------------------------------------------------
  257. c  Computing GAMMA at the cell-centre
  258. c-------------------------------------------------------------
  259.          top=0.0D0
  260.          bot=0.0D0
  261.          do 103  i=1,(nsp-1)
  262.           top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
  263.           bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
  264.  103     continue
  265.          top=cp.gc(nsp)+top
  266.          bot=cv.gc(nsp)+bot
  267.          GAMC=top/bot
  268. C-------------------------------------------------------------
  269.          SEGINI DGDYC
  270.          do 41 i=1,(nsp-1)
  271.            dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  272.      &       GAMC*(cv.gc(i)-cv.gc(nsp)))/bot
  273.  41    continue
  274. c-------------------------------------------------------------
  275. C------------------------------------------------
  276. C******* Densite, vitesse, pression sur le bord
  277. C------------------------------------------------
  278.          ECIN=0.5D0*((UXF*UXF)+(UYF*UYF))
  279.          PSRF=(GM1/GAMF)*(HTF-ECIN)
  280.          RF=PSRF/SF
  281.          RF=RF**(1.0D0/GM1)
  282.          PF=SF*(RF**GAMF)
  283. C------------------------------
  284. C******* Derivatives w.r.t. PC
  285. C------------------------------
  286.          DPSRUN=-1*(GM1/GAMF)*UNC
  287.          DRDUN=USGM1*RF/PSRF*DPSRUN
  288.          DPDUN=GAMF*PSRF*DRDUN
  289.          DUXDUN=CNX
  290.          DUYDUN=CNY
  291. C-------------------------------
  292.          DUNDUX=CNX
  293.          DUNDUY=CNY
  294. C--------------------------------------------------------------
  295. C------------------------------
  296. C******* Derivatives
  297. C------------------------------
  298.         wvec_l(1)=RF
  299.         wvec_l(2)=UXF
  300.         wvec_l(3)=UYF
  301.         wvec_l(4)=PF
  302. C--------------------------
  303.         wvec_r(1)=RC
  304.         wvec_r(2)=UXC
  305.         wvec_r(3)=UYC
  306.         wvec_r(4)=PC
  307. C--------------------------
  308.         nvect(1)=CNX
  309.         nvect(2)=CNY
  310.         tvect(1)=CTX
  311.         tvect(2)=CTY
  312.         call copms2(nsp,jfl,jfr,wvec_l,wvec_r,nvect,tvect,
  313.      &                 yf,mpyc,lrecp,lrecv,nlc)
  314. C--------------------------------------------------------------
  315.         COEF=SURF/VOLU
  316.         JLL=JFL
  317.         JRR=JFR
  318.         SEGACT JLL
  319.         SEGACT JRR
  320.         SEGINI JR
  321. C-------------------------------------------------------------
  322. C   Taking in to account the dependance of the variables
  323. C   at the face on the variables at the centre (UNC)
  324. C-------------------------------------------------------------
  325.        DO 105 I=1,(3+NSP)
  326.        DO 1050 J=1,(3+NSP)
  327.         IF(J .EQ. 2) THEN
  328.           COEF5=(JLL.JAC(I,1)*DRDUN)+(JLL.JAC(I,2)*DUXDUN)+
  329.      &        (JLL.JAC(I,3)*DUYDUN)+(JLL.JAC(I,4)*DPDUN)
  330.           JR.JAC(I,J)=(COEF5*DUNDUX+JRR.JAC(I,J))*COEF
  331.         ELSEIF(J .EQ. 3) THEN
  332.           COEF5=(JLL.JAC(I,1)*DRDUN)+(JLL.JAC(I,2)*DUXDUN)+
  333.      &        (JLL.JAC(I,3)*DUYDUN)+(JLL.JAC(I,4)*DPDUN)
  334.           JR.JAC(I,J)=(COEF5*DUNDUY+JRR.JAC(I,J))*COEF
  335.         ELSE
  336.           JR.JAC(I,J)=JRR.JAC(I,J)*COEF
  337.         ENDIF
  338.  1050  CONTINUE
  339.  105   CONTINUE
  340. C-------------------------------------------------------------
  341. c  matrix wl(i,j) represents the derivative of the i-component
  342. c  of the vector of primitive variables of the left state with
  343. c  respect to the j-component of the vector of the conservative
  344. c  variables of the left state.
  345. c
  346. c  Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  347. c                               vector of primitive variables;
  348. c  (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  349. c                            vector of conservative variables.
  350. c-------------------------------------------------------------
  351.        wl.jac(1,1)=1.0d0
  352.        wl.jac(1,2)=0.0d0
  353.        wl.jac(1,3)=0.0d0
  354.        wl.jac(1,4)=0.0d0
  355.        do 83 i=1,(nsp-1)
  356.          wl.jac(1,4+i)=0.0d0
  357.  83    continue
  358. c------------------------------
  359.        wl.jac(2,1)=-UXC/RC
  360.        wl.jac(2,2)=1.0d0/RC
  361.        wl.jac(2,3)=0.0d0
  362.        wl.jac(2,4)=0.0d0
  363.        do 84 i=1,(nsp-1)
  364.          wl.jac(2,4+i)=0.0d0
  365.  84    continue
  366. c------------------------------
  367.        wl.jac(3,1)=-UYC/RC
  368.        wl.jac(3,2)=0.0d0
  369.        wl.jac(3,3)=1.0d0/RC
  370.        wl.jac(3,4)=0.0d0
  371.        do 85 i=1,(nsp-1)
  372.          wl.jac(3,4+i)=0.0d0
  373.  85    continue
  374. c------------------------------
  375.        br1=0.0d0
  376.        do 86 i=1,(nsp-1)
  377.          br1=br1+dgdyc.vv(i)*yc.yet(i)
  378.  86    continue
  379.        br1=br1*PC/(RC*(GAMC-1.0D0))
  380.        wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1
  381.        wl.jac(4,2)=-UXC*(GAMC-1.0D0)
  382.        wl.jac(4,3)=-UYC*(GAMC-1.0D0)
  383.        wl.jac(4,4)=(GAMC-1.0D0)
  384.        do 87 i=1,(nsp-1)
  385.          wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0))
  386.  87    continue
  387. c------------------------------
  388.        do 88 i=1,(nsp-1)
  389.         do 89 j=1,4
  390.           wl.jac(4+i,j)=0.0d0
  391.           if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC
  392.  89     continue
  393. c------------
  394.         do 890 j=5,(4+nsp-1)
  395.           wl.jac(4+i,j)=0.0d0
  396.           if(4+i.eq.j)  then
  397.             wl.jac(4+i,j)=1.0d0/RC
  398.           endif
  399.  890    continue
  400.  88    continue
  401. c------------------------------------------------
  402. C------------------------------------------------
  403.        do 114 i=1,(3+nsp)
  404.         do 115 j=1,(3+nsp)
  405.          jtt.jac(i,j)=0.0d0
  406.          do 116 k=1,(3+nsp)
  407.           jtt.jac(i,j)=jtt.jac(i,j)+
  408.      &       jr.jac(i,k)*wl.jac(k,j)/COEF
  409.  116     continue
  410.  115    continue
  411.  114  continue
  412. C----------------------------------------------------------------
  413. C******* Jacobian with respect to conservative variables
  414. C----------------------------------------------------------------
  415.          IF(IJAC.EQ.1)THEN
  416.            DO 9 II = 1,(3+NSP)
  417.               DO 15 JJ = 1,(3+NSP)
  418.                 KV = (II-1)*(3+NSP)
  419. C----------------------------------
  420.                 CELL = IMATRI.LIZAFM(1,KV+JJ)
  421.                 CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ)*COEF
  422.  15           CONTINUE
  423.  9         CONTINUE
  424.          ELSEIF(IJAC.EQ.2)THEN
  425.            DO 20 II = 1,(3+NSP)
  426.               DO 25 JJ = 1,(3+NSP)
  427.                 KV = (II-1)*(3+NSP)
  428. C----------------------------------
  429.                 CELL = IMATRI.LIZAFM(1,KV+JJ)
  430.                 CELL.AM(IFAC,1,1) = JR.JAC(II,JJ)
  431.  25           CONTINUE
  432.  20         CONTINUE
  433.          ENDIF
  434. c--------------------------------------------------
  435.       ENDDO
  436. C
  437.       SEGDES MELEFC
  438. C
  439.       SEGDES MLEMC
  440.       SEGDES MLEMCB
  441.       SEGDES MLEMF
  442. C
  443.       SEGDES MPNORM
  444.       SEGDES MPVOL
  445.       SEGDES MPSURF
  446.       SEGDES MPRC
  447.       SEGDES MPPC
  448.       SEGDES MPVC
  449.       SEGDES MPYC
  450.       SEGDES MPLIM
  451.       SEGDES YC
  452.       SEGDES YF
  453.       SEGDES CP
  454.       SEGDES CV
  455.       SEGDES JLL
  456.       SEGDES JRR
  457.       SEGDES JTT
  458.       SEGDES JFL
  459.       SEGDES JFR
  460.       SEGDES JR
  461.       SEGDES WL
  462.       SEGDES  DGDYC
  463.       SEGDES MATRIK
  464.       DO 80 II=1,NBME
  465.        CELL = IMATRI.LIZAFM(1,II)
  466.        SEGDES CELL
  467.  80   CONTINUE
  468.       SEGDES IMATRI
  469. C---------------------------------------------
  470.  9999 CONTINUE
  471.       RETURN
  472.       END
  473.  
  474.  
  475.  
  476.  
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  

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