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

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