Télécharger cli223.eso

Retour à la liste

Numérotation des lignes :

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

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