Télécharger la plaquette Cast3M


fr en

Télécharger cli223.eso

Retour à la liste

Numérotation des lignes : 

cli223

  1. C CLI223    SOURCE    OF166741  24/12/13    21:15:33     12097          
  2.       SUBROUTINE CLI223(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               :  CLI223
  10. C
  11. C DESCRIPTION       :  Subroutine appellée par CLIM22
  12. C                      Jacobian de residu pres de limit
  13. C                      OPTION: 'INRI' 2D
  14. C
  15. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  16. C
  17. C AUTEUR            :  S.Kudriakov, DEN/DM2S/SFME/LTMF
  18. C
  19. C************************************************************************
  20. C
  21. C APPELES (Calcul)  :
  22. C
  23. C************************************************************************
  24. C
  25. C HISTORIQUE (Anomalies et modifications éventuelles)
  26. C
  27. C HISTORIQUE :
  28. C
  29. C************************************************************************
  30. C
  31.       IMPLICIT INTEGER(I-N)
  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,CTX,CTY
  39.      &     ,RF,PF,UXF,UYF
  40.      &     ,UNC,UNF,UTF,SF,ASONC,ASONF,ASON
  41.      &     ,USGM1,G1,G3,ASON2,S,UT,UN,RHO,P,UX,UY
  42.      &     ,DUNDG1,DASDG1,DRHDG1,DPDG1,DUXDG1,DUYDG1
  43.      &     ,DFRDG1,DFMXG1,DFMYG1,DFEDG1
  44.      &     ,DG1DR,DG1DP,DG1DUX,DG1DUY,COEF,DACDG
  45.      &     ,BR1,BOT,TOP,GAMF
  46. c     &     ,EPS,CACCA,CACC1,CEL
  47.  
  48.       CHARACTER*(8) TYPE
  49.  
  50. -INC PPARAM
  51. -INC CCOPTIO
  52. -INC SMLMOTS
  53. -INC SMELEME
  54.       POINTEUR MELEFC.MELEME
  55. -INC SMLENTI
  56.       POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI
  57. -INC SMCHPOI
  58.       POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL,
  59.      &     MPVC.MPOVAL, MPPC.MPOVAL, MPLIM.MPOVAL, MPYC.MPOVAL
  61.       POINTEUR CELL.IZAFM
  62. C-------------------------------------------------------
  63. -INC SMLREEL
  64.       POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  65. C-------------------------------------------------------
  66. C*********  Les Jacobians ******************************
  67. C-------------------------------------------------------
  68.       SEGMENT JACEL
  69.          REAL*8 JAC(3+NSP,3+NSP)
  70.       ENDSEGMENT
  71.       POINTEUR JTL.JACEL, WL.JACEL, JTT.JACEL
  72. C-------------------------------------------------------------
  73. C******* Les fractionines massiques **************************
  74. C-------------------------------------------------------------
  75.        SEGMENT FRAMAS
  76.          REAL*8 YET(NSP)
  77.       ENDSEGMENT
  78.       POINTEUR YC.FRAMAS, YF.FRAMAS
  79. C-------------------------------------------------------
  80. C**********  Les CP's and CV's   ***********************
  81. C-------------------------------------------------------
  82.       SEGMENT GCONST
  83.          REAL*8 GC(NSP)
  84.       ENDSEGMENT
  85.       POINTEUR CP.GCONST, CV.GCONST
  86. C-------------------------------------------------------------
  87. C**********  Segments for the vectors  ***********************
  88. C-------------------------------------------------------------
  89.        SEGMENT VECEL
  90.          REAL*8 VV(NSP)
  91.        ENDSEGMENT
  92.        POINTEUR DYDG1.VECEL, DFRYG1.VECEL,
  93.      &          DG1DY.VECEL, DGDYC.VECEL
  94. C----------------------------------------------------
  95.       SEGINI JTL
  96.       SEGINI JTT
  97.       SEGINI WL
  98. C----------------------------------------------------
  99. C**** KRIPAD pour la correspondance global/local
  100. C----------------------------------------------------
  101.       CALL KRIPAD(MELEMC,MLEMC)
  102.       CALL KRIPAD(MELECB,MLEMCB)
  103.       CALL KRIPAD(MELEMF,MLEMF)
  104. C----------------------------------------------------
  105. C**** CHPOINTs de la table DOMAINE
  106. C----------------------------------------------------
  107.       CALL LICHT(INORM,MPNORM,TYPE,ICEL)
  108.       CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL)
  109.       CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL)
  110. C----------------------------------------------------
  111. C**** CHPOINTs des variables
  112. C----------------------------------------------------
  113.       CALL LICHT(IROC,MPRC,TYPE,ICEL)
  114.       CALL LICHT(IVITC,MPVC,TYPE,ICEL)
  115.       CALL LICHT(IPC,MPPC,TYPE,ICEL)
  116.       CALL LICHT(IYC,MPYC,TYPE,ICEL)
  117.       CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
  118. C--------------------------------------------------------
  119. C**** Boucle sur le face pour le calcul des invariants de
  120. C     Riemann et du flux
  121. C--------------------------------------------------------
  122.       SEGACT MELEFC
  123.       NFAC=MELEFC.NUM(/2)
  124. C---------------------------------
  125. C**** Objet MATRIK
  126. C---------------------------------
  127.       NRIGE = 7
  128.       NMATRI = 1
  129.       NKID = 9
  130.       NKMT = 7
  131. C---------------------------------
  132.       SEGINI MATRIK
  133.       IJACO = MATRIK
  134.       MATRIK.IRIGEL(1,1) = MELRES
  135.       MATRIK.IRIGEL(2,1) = MELRES
  136. C---------------------------------
  137. C**** Matrice non symetrique
  138. C---------------------------------
  139.       MATRIK.IRIGEL(7,1) = 2
  140. C---------------------------------
  141.       NBME = (3+NSP)*(3+NSP)
  142.       NBSOUS = 1
  143.       SEGINI IMATRI
  144.       IF(IJAC.EQ.1)THEN
  145.          MLMOTS=ILIINC
  146.       ELSEIF(IJAC.EQ.2)THEN
  147.          MLMOTS=ILIINP
  148.       ENDIF
  149.       SEGACT MLMOTS
  150.       MATRIK.IRIGEL(4,1) = IMATRI
  151. C-------------------------------------------
  152.       DO 1 J=1,(NSP+3)
  153.         KV=(J-1)*(3+NSP)
  154.         IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1)
  155.         IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2)
  156.         IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3)
  157.         IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4)
  158.         DO 2 I=1,(NSP-1)
  159.           IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I)
  160.  2      CONTINUE
  161.  1    CONTINUE
  162. C-----------------------------------------------
  163.       SEGDES MLMOTS
  164.       MLMOTS=ILIINC
  165.       SEGACT MLMOTS
  166. C-----------------------------------------------
  167.       DO 3 J=1,(NSP+3)
  168.         KV=(J-1)*(3+NSP)
  169.         IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j)
  170.         IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j)
  171.         IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j)
  172.         IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j)
  173.         DO 4 I=1,(NSP-1)
  174.           IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j)
  175.  4     CONTINUE
  176.  3    CONTINUE
  177. C-----------------------------------------------
  178. C-----------------------------------------------
  179.       SEGDES MLMOTS
  180.       NBEL = NFAC
  181.       NBSOUS = 1
  182.       NP = 1
  183.       MP = 1
  184. C-----------------------------------------------------------
  185. C-----------------------------------------------------------
  186.       DO 5 I=1,NBME
  187.         SEGINI CELL
  188.         IMATRI.LIZAFM(1,I) = CELL
  189.  5    CONTINUE
  190. C---------------------------------
  191. C---------------------------------
  192. C---------------------------------
  193. C**** Fin definition MATRIK
  194. C---------------------------------
  195.       DO IFAC=1,NFAC,1
  196.          NGF=MELEFC.NUM(1,IFAC)
  197.          NGC=MELEFC.NUM(2,IFAC)
  198.          NLF=MLEMF.LECT(NGF)
  199.          NLC=MLEMC.LECT(NGC)
  200.          NLCB=MLEMCB.LECT(NGF)
  201.          VOLU=MPVOL.VPOCHA(NLC,1)
  202.          SURF=MPSURF.VPOCHA(NLF,1)
  203. C        In CASTEM les normales sont sortantes
  204.          CNX=-1*MPNORM.VPOCHA(NLF,1)
  205.          CNY=-1*MPNORM.VPOCHA(NLF,2)
  206.          CTX=-1.0D0*CNY
  207.          CTY=CNX
  208. C----------------------------------------------
  209.          SEGINI CP, CV
  210.          MLRECP = LRECP
  211.          MLRECV = LRECV
  212.          SEGACT MLRECP, MLRECV
  213.          DO 10 I=1,(NSP-1)
  214.            CP.GC(I)=MLRECP.PROG(I)
  215.            CV.GC(I)=MLRECV.PROG(I)
  216.  10      CONTINUE
  217.          CP.GC(NSP)=MLRECP.PROG(NSP)
  218.          CV.GC(NSP)=MLRECV.PROG(NSP)
  219. C---------------------------------
  220. C        Variables au centre
  221. C---------------------------------
  222.          RC=MPRC.VPOCHA(NLC,1)
  223.          PC=MPPC.VPOCHA(NLC,1)
  224.          UXC=MPVC.VPOCHA(NLC,1)
  225.          UYC=MPVC.VPOCHA(NLC,2)
  226.          SEGINI YC
  227.          SEGACT MPYC
  228.          DO 100 I=1,(NSP-1)
  229.            YC.YET(I)=MPYC.VPOCHA(NLC,I)
  230.  100     CONTINUE
  231. C---------------------------------
  232. C        Variables à la face
  233. C---------------------------------
  234.          RF=MPLIM.VPOCHA(NLCB,1)
  235.          UXF=MPLIM.VPOCHA(NLCB,2)
  236.          UYF=MPLIM.VPOCHA(NLCB,3)
  237.          PF=MPLIM.VPOCHA(NLCB,IDIM+2)
  238.          SEGINI YF
  239.          SEGACT MPLIM
  240.          DO 101 I=1,(NSP-1)
  241.            YF.YET(I)=MPLIM.VPOCHA(NLCB,IDIM+2+I)
  242.  101     CONTINUE
  243. c-------------------------------------------------------------
  244. c  Computing GAMMA at the cell-center
  245. c-------------------------------------------------------------
  246.          top=0.0D0
  247.          bot=0.0D0
  248.          do 102  i=1,(nsp-1)
  249.           top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
  250.           bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
  251.  102     continue
  252.          top=cp.gc(nsp)+top
  253.          bot=cv.gc(nsp)+bot
  254.          GAMC=top/bot
  255. C-------------------------------------------------------------
  256.          SEGINI DGDYC
  257.          do 41 i=1,(nsp-1)
  258.            dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  259.      &       GAMC*(cv.gc(i)-cv.gc(nsp)))/bot
  260.  41    continue
  261. c-------------------------------------------------------------
  262. c  Computing GAMMA at the face-center
  263. c-------------------------------------------------------------
  264.          top=0.0D0
  265.          bot=0.0D0
  266.          do 103  i=1,(nsp-1)
  267.           top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp))
  268.           bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp))
  269.  103     continue
  270.          top=cp.gc(nsp)+top
  271.          bot=cv.gc(nsp)+bot
  272.          GAMF=top/bot
  273. C--------------------------------------------------------------
  274. C******* On calcule UN, UT ASON, S
  275. C--------------------------------------------------------------
  276.          UNC=(UXC*CNX)+(UYC*CNY)
  277.          UNF=(UXF*CNX)+(UYF*CNY)
  278.          UTF=(UXF*CTX)+(UYF*CTY)
  279. C---------------------------------
  280.          ASONC=(GAMC*PC/RC)**0.5D0
  281.          ASONF=(GAMF*PF/RF)**0.5D0
  282. C---------------------------------
  283.          SF=PF/(RF**GAMF)
  284. C------------------------------------------------
  285. C******* Densite, vitesse, pression sur le bord
  286. C------------------------------------------------
  287.          G1=UNC-(2.0D0*ASONC)/(GAMC-1.0D0)
  288.          G3=UNF+(2.0D0*ASONF)/(GAMF-1.0D0)
  289.          UN=0.5D0*(G1+G3)
  290.          ASON=(0.5D0*(G3-G1))
  291.          ASON=ASON*(GAMF-1.0D0)/2.0D0
  292.          ASON2=ASON*ASON
  293.          S=SF
  294.          UT=UTF
  295.          RHO=ASON2/(GAMF*S)
  296.          RHO=RHO**(1.0D0/(GAMF-1.0D0))
  297.          P=RHO*ASON2/GAMF
  298.          UX=(UN*CNX)+(UT*CTX)
  299.          UY=(UN*CNY)+(UT*CTY)
  300. C------------------------------
  301. C******* Derivatives
  302. C------------------------------
  303.          SEGINI DYDG1
  304.          DUNDG1=0.5D0
  305.          DASDG1=-0.5D0*(GAMF-1.0D0)/2.0D0
  306.          DRHDG1=GAMF*S
  307.          DRHDG1=1.0D0/DRHDG1
  308.          DRHDG1=DRHDG1**(1.0D0/(GAMF-1.0D0))
  309.          DRHDG1=DRHDG1*2.0D0/(GAMF-1.0D0)
  310.          DRHDG1=DRHDG1*(ASON**((3.0D0-GAMF)/(GAMF-1.0D0)))
  311.          DRHDG1=DRHDG1*DASDG1
  312.          DPDG1=((ASON2/GAMF)*DRHDG1)+(((2*ASON*RHO)/GAMF)*DASDG1)
  313.          DO 104  I=1,(NSP-1)
  314.            DYDG1.VV(I) = 0.0D0
  315.  104     CONTINUE
  316. CC------------------------------------------------------
  317. CC******* Test
  318. CC------------------------------------------------------
  319. c         USGM1 = 1.0D0/(GAMF-1.0D0)
  320. c         DFEDG1=(DUNDG1*GAMF*USGM1*P) + (UN*GAMF*USGM1*DPDG1) +
  321. cc     &        (0.5D0*DRHDG1*UN*((UN*UN)+(UT*UT))) +
  322. c     &        (0.5D0*RHO*DUNDG1*((UN*UN)+(UT*UT))) +
  323. c     &        (RHO*UN*UN*DUNDG1)
  324. c         CACCA=UN*GAMF*USGM1*P + 0.5D0*RHO*UN*(UN*UN+UT*UT)
  325. cc---------------
  326. c         DFRYG1.VV(1) = (DRHDG1*UN*YF.YET(1))+
  327. c     &      (RHO*DUNDG1*YF.YET(1))+(RHO*UN*DYDG1.VV(1))
  328. c         CACC1 =RHO*UN*YF.YET(1)
  329. cc---------------
  330. c         EPS=1.0D-8
  331. c         G1=G1*(1+EPS)
  332. c         CEL=UN
  333. c         UN=0.5D0*(G1+G3)
  334. c         write(*,*) DUNDG1, ((UN - CEL)/(EPS*G1))
  335. c         CEL=ASON
  336. c         ASON=(0.5D0*(G3-G1))
  337. c         ASON=ASON*(GAMM-1.0D0)/2.0D0
  338. c         write(*,*) DASDG1, ((ASON - CEL)/(EPS*G1))
  339. c         ASON2=ASON*ASON
  340. c         S=SF
  341. c         UT=UTF
  342. c         CEL=RHO
  343. c         RHO=ASON2/(GAMF*S)
  344. c         RHO=RHO**(1.0D0/(GAMF-1.0D0))
  345. c         write(*,*) DRHDG1, ((RHO - CEL)/(EPS*G1))
  346. c         CEL=P
  347. c         P=RHO*ASON2/GAMF
  348. c         write(*,*) DPDG1, ((P - CEL)/(EPS*G1))
  349. c         CEL=CACCA
  350. c         CACCA=UN*GAMF*USGM1*P + 0.5D0*RHO*UN*(UN*UN+UT*UT)
  351. c         write(*,*) DFEDG1, ((CACCA - CEL)/(EPS*G1))
  352. c         CEL=CACC1
  353. c         CACCA =RHO*UN*YF.YET(1)
  354. c         write(*,*) DFRYG1.VV(1), ((CACCA - CEL)/(EPS*G1))
  355. CC
  356. CC************************************************************
  357. CC
  358.          DUXDG1=DUNDG1*CNX
  359.          DUYDG1=DUNDG1*CNY
  360. C--------------------------------------------------------------
  361. C -------  Be carefull here for outlet BC !!!!!!!!!!!!!!!!!!!!
  362. C--------------------------------------------------------------
  363.          SEGINI DFRYG1
  364.          USGM1 = 1.0D0/(GAMF-1.0D0)
  365.          DFRDG1=(DRHDG1*UN)+(RHO*DUNDG1)
  366.          DFMXG1=(DRHDG1*UN*UX)+(RHO*DUNDG1*UX)+
  367.      &        (RHO*UN*DUXDG1)+(DPDG1*CNX)
  368.          DFMYG1=(DRHDG1*UN*UY)+(RHO*DUNDG1*UY)+
  369.      &        (RHO*UN*DUYDG1)+(DPDG1*CNY)
  370.          DFEDG1=(DUNDG1*GAMF*USGM1*P) + (UN*GAMF*USGM1*DPDG1) +
  371.      &        (0.5D0*DRHDG1*UN*((UN*UN)+(UT*UT))) +
  372.      &        (0.5D0*RHO*DUNDG1*((UN*UN)+(UT*UT))) +
  373.      &        (RHO*UN*UN*DUNDG1)
  374.          DO 105 I=1,(NSP-1)
  375.            DFRYG1.VV(I) = (DRHDG1*UN*YF.YET(I))+
  376.      &      (RHO*DUNDG1*YF.YET(I))+(RHO*UN*DYDG1.VV(I))
  377.  105     CONTINUE
  378. C-----------------------------------------------------
  379. C******* Jacobian with respect to primitive variables
  380. C-----------------------------------------------------
  381.          SEGINI DG1DY
  382.          USGM1 = 1.0D0/(GAMC-1.0D0)
  383.          DG1DR=USGM1*ASONC/RC
  384.          DG1DP=-1.0D0*USGM1*ASONC/PC
  385.          DG1DUX=CNX
  386.          DG1DUY=CNY
  387.          DACDG=0.5D0*ASONC/GAMC
  388.          DO 106 I=1,(NSP-1)
  389.           DG1DY.VV(I)=2.0D0*USGM1*(ASONC*USGM1-DACDG)*
  390.      &               dgdyc.vv(i)
  391.  106     CONTINUE
  392. C----------------------------------------
  393.          COEF=SURF/VOLU
  394. C----------------------------------------
  395.          JTL.JAC(1,1) = DFRDG1*DG1DR*COEF
  396.          JTL.JAC(1,2) = DFRDG1*DG1DUX*COEF
  397.          JTL.JAC(1,3) = DFRDG1*DG1DUY*COEF
  398.          JTL.JAC(1,4) = DFRDG1*DG1DP*COEF
  399.          DO 107 I=1,(NSP-1)
  400.            JTL.JAC(1,4+I) = DFRDG1*DG1DY.VV(I)*COEF
  401.  107     CONTINUE
  402. c         CACCA = G1
  403. c          top=0.0D0
  404. c          bot=0.0D0
  405. c          yc.yet(1) = yc.yet(1)+EPS
  406. c         do 1023  i=1,(nsp-1)
  407. c          top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
  408. c          bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
  409. c 1023  continue
  410. c          top=cp.gc(nsp)+top
  411. c          bot=cv.gc(nsp)+bot
  412. c          GAMC=top/bot
  413. c          ASONC=(GAMC*PC/RC)**(0.5D0)
  414. c          GAMM = 0.5D0*(GAMC+GAMF)
  415. c          G1=UNC-(2.0D0*ASONC)/(GAMC-1.0D0)
  416. c          G3=UNF+(2.0D0*ASONF)/(GAMF-1.0D0)
  417. c          UN=0.5D0*(G1+G3)
  418. c          ASON=(0.5D0*(G3-G1))
  419. c          ASON=ASON*(GAMM-1.0D0)/2.0D0
  420. c          ASON2=ASON*ASON
  421. c          S=SF
  422. c          UT=UTF
  423. c          RHO=ASON2/(GAMF*S)
  424. c          RHO=RHO**(1.0D0/(GAMF-1.0D0))
  425. c          CEL=CACCA
  426. c          CACCA=G1
  427. C----------------------------------------
  428.          JTL.JAC(2,1) = DFMXG1*DG1DR*COEF
  429.          JTL.JAC(2,2) = DFMXG1*DG1DUX*COEF
  430.          JTL.JAC(2,3) = DFMXG1*DG1DUY*COEF
  431.          JTL.JAC(2,4) = DFMXG1*DG1DP*COEF
  432.          DO 108 I=1,(NSP-1)
  433.            JTL.JAC(2,4+I) = DFMXG1*DG1DY.VV(I)*COEF
  434.  108     CONTINUE
  435. C----------------------------------------
  436.          JTL.JAC(3,1) = DFMYG1*DG1DR*COEF
  437.          JTL.JAC(3,2) = DFMYG1*DG1DUX*COEF
  438.          JTL.JAC(3,3) = DFMYG1*DG1DUY*COEF
  439.          JTL.JAC(3,4) = DFMYG1*DG1DP*COEF
  440.          DO 109 I=1,(NSP-1)
  441.            JTL.JAC(3,4+I) = DFMYG1*DG1DY.VV(I)*COEF
  442.  109     CONTINUE
  443. C----------------------------------------
  444.          JTL.JAC(4,1) = DFEDG1*DG1DR*COEF
  445.          JTL.JAC(4,2) = DFEDG1*DG1DUX*COEF
  446.          JTL.JAC(4,3) = DFEDG1*DG1DUY*COEF
  447.          JTL.JAC(4,4) = DFEDG1*DG1DP*COEF
  448.          DO 110 I=1,(NSP-1)
  449.            JTL.JAC(4,4+I) = DFEDG1*DG1DY.VV(I)*COEF
  450.  110     CONTINUE
  451. C----------------------------------------
  452.          DO 111 I=1,(NSP-1)
  453.           JTL.JAC(4+I,1) = DFRYG1.VV(I)*DG1DR*COEF
  454.           JTL.JAC(4+I,2) = DFRYG1.VV(I)*DG1DUX*COEF
  455.           JTL.JAC(4+I,3) = DFRYG1.VV(I)*DG1DUY*COEF
  456.           JTL.JAC(4+I,4) = DFRYG1.VV(I)*DG1DP*COEF
  457.           DO 112 J=1,(NSP-1)
  458.            JTL.JAC(4+I,4+J) = DFRYG1.VV(I)*DG1DY.VV(J)*COEF
  459.  112      CONTINUE
  460.  111     CONTINUE
  461. C---------------------------------------------------------
  462. c  matrix wl(i,j) represents the derivative of the i-component
  463. c  of the vector of primitive variables of the left state with
  464. c  respect to the j-component of the vector of the conservative
  465. c  variables of the left state.
  466. c
  467. c  Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  468. c                               vector of primitive variables;
  469. c  (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  470. c                            vector of conservative variables.
  471. c-------------------------------------------------------------
  472.        wl.jac(1,1)=1.0d0
  473.        wl.jac(1,2)=0.0d0
  474.        wl.jac(1,3)=0.0d0
  475.        wl.jac(1,4)=0.0d0
  476.        do 83 i=1,(nsp-1)
  477.          wl.jac(1,4+i)=0.0d0
  478.  83    continue
  479. c------------------------------
  480.        wl.jac(2,1)=-UXC/RC
  481.        wl.jac(2,2)=1.0d0/RC
  482.        wl.jac(2,3)=0.0d0
  483.        wl.jac(2,4)=0.0d0
  484.        do 84 i=1,(nsp-1)
  485.          wl.jac(2,4+i)=0.0d0
  486.  84    continue
  487. c------------------------------
  488.        wl.jac(3,1)=-UYC/RC
  489.        wl.jac(3,2)=0.0d0
  490.        wl.jac(3,3)=1.0d0/RC
  491.        wl.jac(3,4)=0.0d0
  492.        do 85 i=1,(nsp-1)
  493.          wl.jac(3,4+i)=0.0d0
  494.  85    continue
  495. c------------------------------
  496.        br1=0.0d0
  497.        do 86 i=1,(nsp-1)
  498.          br1=br1+dgdyc.vv(i)*yc.yet(i)
  499.  86    continue
  500.        br1=br1*PC/(RC*(GAMC-1.0D0))
  501.        wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1
  502.        wl.jac(4,2)=-UXC*(GAMC-1.0D0)
  503.        wl.jac(4,3)=-UYC*(GAMC-1.0D0)
  504.        wl.jac(4,4)=(GAMC-1.0D0)
  505.        do 87 i=1,(nsp-1)
  506.          wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0))
  507.  87    continue
  508. c------------------------------
  509.        do 88 i=1,(nsp-1)
  510.         do 89 j=1,4
  511.           wl.jac(4+i,j)=0.0d0
  512.           if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC
  513.  89     continue
  514. c------------
  515.         do 890 j=5,(4+nsp-1)
  516.           wl.jac(4+i,j)=0.0d0
  517.           if(4+i.eq.j)  then
  518.             wl.jac(4+i,j)=1.0d0/RC
  519.           endif
  520.  890    continue
  521.  88    continue
  522. c------------------------------------------------
  523. C------------------------------------------------
  524.        do 114 i=1,(3+nsp)
  525.         do 115 j=1,(3+nsp)
  526.          jtt.jac(i,j)=0.0d0
  527.          do 116 k=1,(3+nsp)
  528.           jtt.jac(i,j)=jtt.jac(i,j)+jtl.jac(i,k)*wl.jac(k,j)
  529.  116     continue
  530.  115    continue
  531.  114  continue
  532. C----------------------------------------------------------------
  533. C******* Jacobian with respect to conservative variables
  534. C----------------------------------------------------------------
  535.          IF(IJAC.EQ.1)THEN
  536.            DO 9 II = 1,(3+NSP)
  537.               DO 15 JJ = 1,(3+NSP)
  538.                 KV = (II-1)*(3+NSP)
  539. C----------------------------------
  540.                 CELL = IMATRI.LIZAFM(1,KV+JJ)
  541.                 CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ)
  542.  15           CONTINUE
  543.  9         CONTINUE
  544.          ELSEIF(IJAC.EQ.2)THEN
  545.            DO 20 II = 1,(3+NSP)
  546.               DO 25 JJ = 1,(3+NSP)
  547.                 KV = (II-1)*(3+NSP)
  548. C----------------------------------
  549.                 CELL = IMATRI.LIZAFM(1,KV+JJ)
  550.                 CELL.AM(IFAC,1,1) = JTL.JAC(II,JJ)
  551.  25           CONTINUE
  552.  20         CONTINUE
  553.          ENDIF
  554. c--------------------------------------------------
  555.       ENDDO
  556. C
  557.       SEGDES MELEFC
  558. C
  559.       SEGSUP MLEMC
  560.       SEGSUP MLEMCB
  561.       SEGSUP MLEMF
  562. C
  563.       SEGDES MPNORM
  564.       SEGDES MPVOL
  565.       SEGDES MPSURF
  566.       SEGDES MPRC
  567.       SEGDES MPPC
  568.       SEGDES MPVC
  569.       SEGDES MPYC
  570.       SEGDES MPLIM
  571.       SEGDES YC
  572.       SEGDES YF
  573.       SEGDES CP
  574.       SEGDES CV
  575.       SEGDES JTL
  576.       SEGDES JTT
  577.       SEGDES WL
  578.       SEGDES  DYDG1, DFRYG1,
  579.      &          DG1DY, DGDYC
  580.       SEGDES MATRIK
  581.       DO 80 II=1,NBME
  582.        CELL = IMATRI.LIZAFM(1,II)
  583.        SEGDES CELL
  584.  80   CONTINUE
  585.       SEGDES IMATRI
  586. C---------------------------------------------
  587.  9999 CONTINUE
  588.       RETURN
  589.       END
  590.  
  591.  
  592.  
  593.  
  594.  
  595.  
  596.  
  597.  
  598.  
  599.  
  600.  
  601.  
  602.  
  603.  
  604.  
  605.  

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