Télécharger konja5.eso

Retour à la liste

Numérotation des lignes :

  1. C KONJA5 SOURCE PV 16/11/17 21:59:57 9180
  2. SUBROUTINE KONJA5(IIMPL,ILINC,IROF,IPF,IVITF,IGAMF,
  3. & IRN,IPN,IVN,IGAMN,
  4. $ ICHPVO,ICHPSU, ICHPNO, MELEMC, MELEFE, MELEMF, MELTFA,
  5. & IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN,
  6. & ICRN, ICVN, ICPN, IVLIM,MELLIM,IMAT)
  7. C
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : KONJA5
  13. C
  14. C DESCRIPTION : Voir KONV11
  15. C Calcul du jacobien du résidu "2nd order accurate"
  16. C
  17. C Cas deux dimensions, gaz "calorically perfect"
  18. C
  19. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  20. C
  21. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF
  22. C
  23. C************************************************************************
  24. C
  25. C APPELES (Outils
  26. C CASTEM) : KRIPAD, LICHT, ERREUR
  27. C
  28. C APPELES (Calcul) : VLHJ3, AUSMW, VLHJ1, AUSMF
  29. C
  30. C************************************************************************
  31. C
  32. C ENTREES
  33. C
  34. C IIMPL : INTEGER
  35. C = 4 -> methode VLH
  36. C = 5 -> methode AUSM+
  37. C
  38. C ILINC : liste des inconnues (pointeur d'un objet de type
  39. C LISTMOTS)
  40. C
  41. C 1) Pointeurs des CHPOINT/CHAMELEM
  42. C
  43. C IROF : MCHAML 'FACEL' densité
  44. C
  45. C IPF : MCHAML 'FACEL' pression
  46. C
  47. C IVITF : MCHAML 'FACEL' vitesse/normales
  48. C
  49. C IGAMF : MCHAML 'FACEL' gamma
  50. C
  51. C IRN : CHPOINT 'CENTRE' densité
  52. C
  53. C IPN : CHPOINT 'CENTRE' pression
  54. C
  55. C IVN : CHPOINT 'CENTRE' vitesse
  56. C
  57. C IGAMN : CHPOINT 'CENTRE' gamma
  58. C
  59. C ICHPVO : CHPOINT VOLUME contenant le volume
  60. C
  61. C ICHPSU : CHPOINT FACE contenant la surface des faces
  62. C
  63. C ICHPNO : CHPOINT 'FACE' contenant les normales aux faces
  64. C
  65. C IGRN : CHPOINT 'CENTRE' gradient de densité
  66. C
  67. C IGVN : CHPOINT 'CENTRE' gradient de vitesse
  68. C
  69. C IGPN : CHPOINT 'CENTRE' gradient de pression
  70. C
  71. C IGLRN : CHPOINT 'CENTRE' limiteur du gradient de densité
  72. C
  73. C IGLVN : CHPOINT 'CENTRE' limiteur du gradient de vitesse
  74. C
  75. C IGLPN : CHPOINT 'CENTRE' limiteur du gradient de pression
  76. C
  77. C ICRN : MCHAML contenant les coeff. pour calculer le gradient
  78. C de densité
  79. C
  80. C ICVN : MCHAML contenant les coeff. pour calculer le gradient
  81. C de vitesse
  82. C
  83. C ICPN : MCHAML contenant les coeff. pour calculer le gradient
  84. C de pression
  85. C
  86. C IVLIM : CHPOINT des conditions aux limites sur la vitesse
  87. C
  88. C 2) Pointeurs de MELEME de la table DOMAINE
  89. C
  90. C MELEMC : MELEME 'CENTRE', SPG des CENTRES des elts
  91. C
  92. C MELEFE : MELEME 'FACEL' (centre elt gauche, centre de face
  93. C centre elt droite)
  94. C
  95. C MELEMF : MELEME 'FACE', SPG des CENTRES des faces
  96. C
  97. C MELTFA : MELEME 'ELTFA' (centre d'elt, centres de faces)
  98. C
  99. C MELLIM : MELEME SPG des conditions aux bords
  100. C
  101. C SORTIES
  102. C
  103. C IMAT : pointeur de la MATRIK du jacobien du residu
  104. C
  105. C************************************************************************
  106. C
  107. C HISTORIQUE (Anomalies et modifications éventuelles)
  108. C
  109. C HISTORIQUE :
  110. C
  111. C************************************************************************
  112. C
  113. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  114. C GAMMA \in (1,3)
  115. C Si non il faut le faire!!!
  116. C
  117. C************************************************************************
  118. C
  119. C
  120. C**** Variables de COOPTIO
  121. C
  122. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  123. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  124. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  125. C & ,IECHO, IIMPI, IOSPI
  126. C & ,IDIM, IFICLE, IPREFI
  127. CC & ,MCOORD
  128. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  129. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  130. C & ,NORINC,NORVAL,NORIND,NORVAD
  131. C & ,NUCROU, IPSAUV
  132. C
  133. IMPLICIT INTEGER(I-N)
  134. INTEGER IIMPL,ILINC, ICHPVO, ICHPNO, ICHPSU
  135. & ,IROF,IVITF,IPF,IGAMF, IRN,IPN,IVN,IGAMN
  136. & , IGRN, IGVN, IGPN,IGLRN, IGLVN, IGLPN
  137. & , ICRN, ICVN, ICPN, IVLIM , IMAT
  138. & , NBSOUP, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID
  139. & , NKMT, NBME, NBEL, MP, NP, NP1
  140. & , IFAC, NGCF, NLCF, NGCF1, NGCG, NGCD
  141. & , NGC, NLC, NLCD, NGNO, NGN1,NGN2, NLNO, IELEM, IELPRI
  142. & , JG, NBSOUS, NBSOU1, ISOUS, NBEL1
  143. & , ISPG1, ISPG2
  144. & , I1, IGEOMC, IGEOMF, IGEOML
  145. & , IFACGR, NGGRA, NLGRA, NGFGR, NLFGR, NG, ND, ICOOR, IFAC2
  146. & , NLVIMP
  147. & , NGPRI,NLPRI,IELETR,IDUATR,IPRITR, NLFL
  148. REAL*8 ROG, PG, UNG, UTG, UXG, UYG, RETG, GAMG, VOLG, VOLD
  149. & , ROD, PD, UND, UTD, UXD, UYD
  150. & , SURF, CNX, CNY, CTX, CTY, ECIN, FUNCEL
  151. & , DFRO(4), DFRET(4), DFRUN(4), DFRUT(4),DFRUX(4),DFRUY(4)
  152. & , ALR, ALP, ALVX, ALVY, XCG, YCG, DXCF, DYCF, CELLR, CELLP
  153. & , CVXVX, CVXVY, CVYVX, CVYVY, CNX1, CNY1, CTX1, CTY1
  154. & , DDP, DDRO, DDUX, DDUY
  155. & , GAMGM1, CELCO1, CELCO2, CELCO3, CELCO4, CELCO5
  156.  
  157. CHARACTER*8 TYPE
  158. C
  159. C**** LES INCLUDES
  160. C
  161. -INC SMCOORD
  162. -INC CCOPTIO
  163. -INC SMCHAML
  164. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
  165. & MELT1X.MELVAL, MELT1Y.MELVAL
  166. POINTEUR MELVUN.MELVAL, MELVUT.MELVAL
  167. POINTEUR MELRO.MELVAL, MELP.MELVAL,
  168. & MELGAM.MELVAL
  169. POINTEUR MELCR1.MELVAL, MELCR2.MELVAL, MELCP1.MELVAL,
  170. & MELCP2.MELVAL, MELCV1.MELVAL, MELCV2.MELVAL
  171. C
  172. -INC SMCHPOI
  173. -INC SMELEME
  174. -INC SMLMOTS
  175. -INC SMLENTI
  176. POINTEUR MPVOLU.MPOVAL, MPNORM.MPOVAL, MPOVSU.MPOVAL,
  177. $ MPVGRN.MPOVAL,MPGLRN.MPOVAL, MPVGPN.MPOVAL, MPGLPN.MPOVAL,
  178. $ MPVGVN.MPOVAL,MPGLVN.MPOVAL, MPVLIM.MPOVAL,MPOVRN.MPOVAL,MPOVVN
  179. $ .MPOVAL,MPOVPN.MPOVAL,MPOGAM.MPOVAL
  180. POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME,
  181. & MELTFA.MELEME,MELPRI.MELEME,MELGRR.MELEME,MELGRP.MELEME
  182. $ ,MELGRV.MELEME, MELLIM.MELEME
  183. POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLESPG.MLENTI,MLEPRI.MLENTI
  184. & ,MLEVL.MLENTI,MLEGRR.MLENTI,MLEGRP.MLENTI,MLEGRV.MLENTI
  185. & ,MLELIM.MLENTI
  186. POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RRET.IZAFM,
  187. & UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXRET.IZAFM,
  188. & UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYRET.IZAFM,
  189. & RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETRET.IZAFM
  190. POINTEUR MLMINC.MLMOTS
  191. C
  192. C**** KRIPAD pour la correspondance global/local des conditions limits
  193. C
  194. CALL KRIPAD(MELLIM,MLELIM)
  195. C SEGACT MELLIM
  196. C
  197. C**** Masse volumique
  198. C
  199. MCHEL1 = IROF
  200. SEGACT MCHEL1
  201. MCHAM1 = MCHEL1.ICHAML(1)
  202. SEGACT MCHAM1
  203. MELRO = MCHAM1.IELVAL(1)
  204. SEGDES MCHEL1
  205. SEGDES MCHAM1
  206. C
  207. C**** Pression
  208. C
  209. MCHEL1 = IPF
  210. SEGACT MCHEL1
  211. MCHAM1 = MCHEL1.ICHAML(1)
  212. SEGACT MCHAM1
  213. MELP = MCHAM1.IELVAL(1)
  214. SEGDES MCHEL1
  215. SEGDES MCHAM1
  216. C
  217. C**** Gamma
  218. C
  219. MCHEL1 = IGAMF
  220. SEGACT MCHEL1
  221. MCHAM1 = MCHEL1.ICHAML(1)
  222. SEGACT MCHAM1
  223. MELGAM = MCHAM1.IELVAL(1)
  224. SEGDES MCHEL1
  225. SEGDES MCHAM1
  226. C
  227. C**** Vitesse et cosinus directeurs du repere (n,t)
  228. C
  229. MCHEL1 = IVITF
  230. SEGACT MCHEL1
  231. C
  232. C**** La vitesse a comme SPG MELEFE
  233. C Le cosinus directeurs ont comme SPG MELEMF
  234. C
  235. C MCHAM1 -> Cosinus directeurs
  236. C MCHAM2 -> Vitesse
  237. C
  238. ISPG1 = MCHEL1.IMACHE(1)
  239. ISPG2 = MCHEL1.IMACHE(2)
  240. IF((ISPG1 .EQ. MELEMF) .AND. (ISPG2 .EQ. MELEFE))THEN
  241. MCHAM1 = MCHEL1.ICHAML(1)
  242. MCHAM2 = MCHEL1.ICHAML(2)
  243. ELSEIF((ISPG1 .EQ. MELEFE) .AND. (ISPG2 .EQ. MELEMF))THEN
  244. MCHAM1 = MCHEL1.ICHAML(2)
  245. MCHAM2 = MCHEL1.ICHAML(1)
  246. ELSE
  247. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  248. CALL ERREUR(5)
  249. GOTO 9999
  250. ENDIF
  251. SEGACT MCHAM1
  252. MELVNX = MCHAM1.IELVAL(1)
  253. MELVNY = MCHAM1.IELVAL(2)
  254. MELT1X = MCHAM1.IELVAL(3)
  255. MELT1Y = MCHAM1.IELVAL(4)
  256. SEGDES MCHAM1
  257. SEGACT MCHAM2
  258. MELVUN = MCHAM2.IELVAL(1)
  259. MELVUT = MCHAM2.IELVAL(2)
  260. SEGDES MCHAM2
  261. SEGDES MCHEL1
  262. C
  263. C**** Activation des MCHAMLs
  264. C
  265. SEGACT MELRO
  266. SEGACT MELP
  267. SEGACT MELGAM
  268. SEGACT MELVUN
  269. SEGACT MELVUT
  270. SEGACT MELVNX
  271. SEGACT MELVNY
  272. SEGACT MELT1X
  273. SEGACT MELT1Y
  274. C
  275. C**************** CHPOINT gradients, limiteurs, vitesse aux bords
  276. C
  277. CALL LICHT(IRN,MPOVRN,TYPE,IGEOMC)
  278. IF(IERR.NE.0)GOTO 9999
  279. C SEGACT MPOVRN
  280. CALL LICHT(IVN,MPOVVN,TYPE,IGEOMC)
  281. IF(IERR.NE.0)GOTO 9999
  282. C SEGACT MPOVVN
  283. CALL LICHT(IPN,MPOVPN,TYPE,IGEOMC)
  284. IF(IERR.NE.0)GOTO 9999
  285. C SEGACT MPOVPN
  286. CALL LICHT(IGAMN,MPOGAM,TYPE,IGEOMC)
  287. IF(IERR.NE.0)GOTO 9999
  288. C SEGACT MPOGAM
  289. CALL LICHT(IGRN,MPVGRN,TYPE,IGEOMC)
  290. IF(IERR.NE.0)GOTO 9999
  291. C SEGACT MPVGRN
  292. CALL LICHT(IGLRN,MPGLRN,TYPE,IGEOMC)
  293. IF(IERR.NE.0)GOTO 9999
  294. C SEGACT MPGLRN
  295. CALL LICHT(IGPN,MPVGPN,TYPE,IGEOMC)
  296. IF(IERR.NE.0)GOTO 9999
  297. C SEGACT MPVGPN
  298. CALL LICHT(IGLPN,MPGLPN,TYPE,IGEOMC)
  299. IF(IERR.NE.0)GOTO 9999
  300. C SEGACT MPGLPN
  301. CALL LICHT(IGVN,MPVGVN,TYPE,IGEOMC)
  302. IF(IERR.NE.0)GOTO 9999
  303. C SEGACT MPVGPN
  304. CALL LICHT(IGLVN,MPGLVN,TYPE,IGEOMC)
  305. IF(IERR.NE.0)GOTO 9999
  306. C SEGACT MPGLVN
  307. CALL LICHT(IVLIM,MPVLIM,TYPE,IGEOML)
  308. IF(IERR.NE.0)GOTO 9999
  309. C SEGACT MPVLIM
  310. C
  311. C**************** KRIPAD pour la correspondance global/local des centres d'elts
  312. C
  313. CALL KRIPAD(MELEMC,MLENTC)
  314. IF(IERR.NE.0)GOTO 9999
  315. C
  316. C SEGACT MLENTC
  317. SEGACT MELEMC
  318. C
  319. C**** KRIPAD pour la correspondance global/local des centres des faces
  320. C
  321. CALL KRIPAD(MELEMF,MLENTF)
  322. IF(IERR.NE.0)GOTO 9999
  323. C
  324. C SEGACT MLENTF
  325. C
  326. C**** KRIPAD pour la correspondance global/local des B.C. sur le vitesse
  327. C
  328. CALL KRIPAD(IGEOML,MLEVL)
  329. IF(IERR.NE.0)GOTO 9999
  330. C
  331. C SEGACT MLEVL
  332. C
  333. SEGACT MELEFE
  334. C
  335. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  336. IF(IERR.NE.0)GOTO 9999
  337. CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC)
  338. IF(IERR.NE.0)GOTO 9999
  339. CALL LICHT(ICHPNO,MPNORM,TYPE,IGEOMC)
  340. IF(IERR.NE.0)GOTO 9999
  341. C
  342. C**** LICHT active les MPOVALs en *MOD
  343. C
  344. C i.e.
  345. C
  346. C SEGACT MPOVSU*MOD
  347. C SEGACT MPVOLU*MOD
  348. C SEGACT MPNORM*MOD
  349. C
  350. C
  351. C**** Les gradients ont tous le meme SPG (mais le numero
  352. C du pointeur peut-etre different).
  353. C Ils ont le meme nombre des SPGs que MELTFA
  354. C En plus le SPGs ont le meme ordre:
  355. C le iem elt de MELTFA correspond a l'elt qui a
  356. C comme centre le iem noeud de CENTRE
  357. C
  358. SEGACT MELTFA
  359. NBSOUP=MELTFA.LISOUS(/1)
  360. IF(NBSOUP.EQ.0) NBSOUP=1
  361. JG=NBSOUP
  362. SEGINI MLESPG
  363. IF(NBSOUP .EQ. 1)THEN
  364. MLESPG.LECT(1)=MELTFA
  365. ELSE
  366. DO I1 = 1, NBSOUP, 1
  367. MLESPG.LECT(I1)=MELTFA.LISOUS(I1)
  368. ENDDO
  369. ENDIF
  370. C
  371. MCHEL1=ICRN
  372. SEGACT MCHEL1
  373. NBSOU1=MCHEL1.IMACHE(/1)
  374. IF(NBSOU1 .NE. NBSOUP)THEN
  375. WRITE(IOIMP,*) 'SPG coeffs grads = ???'
  376. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  377. CALL ERREUR(5)
  378. GOTO 9999
  379. ENDIF
  380. C
  381. MCHEL2=ICPN
  382. SEGACT MCHEL2
  383. NBSOU1=MCHEL2.IMACHE(/1)
  384. IF(NBSOU1 .NE. NBSOUP)THEN
  385. WRITE(IOIMP,*) 'SPG coeffs grads = ???'
  386. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  387. CALL ERREUR(5)
  388. GOTO 9999
  389. ENDIF
  390. C
  391. MCHEL2=ICVN
  392. SEGACT MCHEL2
  393. NBSOU1=MCHEL2.IMACHE(/1)
  394. IF(NBSOU1 .NE. NBSOUP)THEN
  395. WRITE(IOIMP,*) 'SPG coeffs grads = ???'
  396. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  397. CALL ERREUR(5)
  398. GOTO 9999
  399. ENDIF
  400. C
  401. C**** N.B.: ICRN, ICPN, ICVN sont activés!!!
  402. C
  403. C
  404. C****************** L'objet matrik
  405. C
  406. C
  407. NRIGE = 7
  408. NMATRI = 1
  409. NKID = 9
  410. NKMT = 7
  411. C
  412. SEGINI MATRIK
  413. IMAT = MATRIK
  414. C
  415. C**** Maillage des inconnues duales: MELEMC
  416. C Maillage des inconnues primales:
  417. C même structure que MCHEL1.IMACHE(ISOUS)
  418. C Donc on partitionne les deux dans la meme
  419. C maniere.
  420. C
  421. JG=NBSOUP
  422. SEGINI MLEPRI
  423. SEGINI MLEGRR
  424. SEGINI MLEGRP
  425. SEGINI MLEGRV
  426. IF(NBSOUP .EQ. 1)THEN
  427. MCHEL1=ICRN
  428. IPT1=MCHEL1.IMACHE(1)
  429. SEGACT IPT1
  430. SEGINI, MELPRI = IPT1
  431. MLEPRI.LECT(1) = MELPRI
  432. MLEGRR.LECT(1) = IPT1
  433. C
  434. MCHEL1=ICPN
  435. IPT1=MCHEL1.IMACHE(1)
  436. SEGACT IPT1
  437. MLEGRP.LECT(1) = IPT1
  438. C
  439. MCHEL1=ICVN
  440. IPT1=MCHEL1.IMACHE(1)
  441. SEGACT IPT1
  442. MLEGRV.LECT(1) = IPT1
  443. ELSE
  444. NBSOUS=NBSOUP
  445. NBREF=0
  446. NBNN=0
  447. NBELEM=0
  448. SEGINI MELPRI
  449. MCHEL1=ICRN
  450. MCHEL2=ICPN
  451. MCHEL3=ICVN
  452. DO ISOUS = 1, NBSOUP, 1
  453. IPT1=MCHEL1.IMACHE(ISOUS)
  454. IPT2=MCHEL2.IMACHE(ISOUS)
  455. IPT3=MCHEL3.IMACHE(ISOUS)
  456. SEGACT IPT1
  457. SEGACT IPT2
  458. SEGACT IPT3
  459. SEGINI, MELEME=IPT1
  460. MELPRI.LISOUS(ISOUS) = MELEME
  461. MLEPRI.LECT(ISOUS) = MELEME
  462. MLEGRR.LECT(ISOUS) = IPT1
  463. MLEGRP.LECT(ISOUS) = IPT2
  464. MLEGRV.LECT(ISOUS) = IPT3
  465. ENDDO
  466. ENDIF
  467. MATRIK.IRIGEL(2,1) = MELPRI
  468. MATRIK.IRIGEL(1,1) = MELPRI
  469. C
  470. C**** Matrice non symetrique
  471. C
  472. MATRIK.IRIGEL(7,1) = 2
  473. NBME = 16
  474. NBSOUS=NBSOUP
  475. SEGINI IMATRI
  476. MLMINC = ILINC
  477. SEGACT MLMINC
  478. MATRIK.IRIGEL(4,1) = IMATRI
  479. C
  480. IMATRI.LISPRI(1) = MLMINC.MOTS(1)
  481. IMATRI.LISPRI(2) = MLMINC.MOTS(2)
  482. IMATRI.LISPRI(3) = MLMINC.MOTS(3)
  483. IMATRI.LISPRI(4) = MLMINC.MOTS(4)
  484. IMATRI.LISPRI(5) = MLMINC.MOTS(1)
  485. IMATRI.LISPRI(6) = MLMINC.MOTS(2)
  486. IMATRI.LISPRI(7) = MLMINC.MOTS(3)
  487. IMATRI.LISPRI(8) = MLMINC.MOTS(4)
  488. IMATRI.LISPRI(9) = MLMINC.MOTS(1)
  489. IMATRI.LISPRI(10) = MLMINC.MOTS(2)
  490. IMATRI.LISPRI(11) = MLMINC.MOTS(3)
  491. IMATRI.LISPRI(12) = MLMINC.MOTS(4)
  492. IMATRI.LISPRI(13) = MLMINC.MOTS(1)
  493. IMATRI.LISPRI(14) = MLMINC.MOTS(2)
  494. IMATRI.LISPRI(15) = MLMINC.MOTS(3)
  495. IMATRI.LISPRI(16) = MLMINC.MOTS(4)
  496. C
  497. IMATRI.LISDUA(1) = MLMINC.MOTS(1)
  498. IMATRI.LISDUA(2) = MLMINC.MOTS(1)
  499. IMATRI.LISDUA(3) = MLMINC.MOTS(1)
  500. IMATRI.LISDUA(4) = MLMINC.MOTS(1)
  501. IMATRI.LISDUA(5) = MLMINC.MOTS(2)
  502. IMATRI.LISDUA(6) = MLMINC.MOTS(2)
  503. IMATRI.LISDUA(7) = MLMINC.MOTS(2)
  504. IMATRI.LISDUA(8) = MLMINC.MOTS(2)
  505. IMATRI.LISDUA(9) = MLMINC.MOTS(3)
  506. IMATRI.LISDUA(10) = MLMINC.MOTS(3)
  507. IMATRI.LISDUA(11) = MLMINC.MOTS(3)
  508. IMATRI.LISDUA(12) = MLMINC.MOTS(3)
  509. IMATRI.LISDUA(13) = MLMINC.MOTS(4)
  510. IMATRI.LISDUA(14) = MLMINC.MOTS(4)
  511. IMATRI.LISDUA(15) = MLMINC.MOTS(4)
  512. IMATRI.LISDUA(16) = MLMINC.MOTS(4)
  513. C
  514. C**** Boucle sur les LISOUS de MELPRI
  515. C
  516. IELEM=0
  517. C
  518. DO ISOUS=1,NBSOUP,1
  519. C
  520. C******* Les coeffs des grad.
  521. C
  522. MCHEL1=ICRN
  523. MCHAM1=MCHEL1.ICHAML(ISOUS)
  524. SEGACT MCHAM1
  525. MELCR1=MCHAM1.IELVAL(1)
  526. MELCR2=MCHAM1.IELVAL(2)
  527. SEGACT MELCR1, MELCR2
  528. SEGDES MCHAM1
  529. C
  530. MCHEL1=ICPN
  531. MCHAM1=MCHEL1.ICHAML(ISOUS)
  532. SEGACT MCHAM1
  533. MELCP1=MCHAM1.IELVAL(1)
  534. MELCP2=MCHAM1.IELVAL(2)
  535. SEGACT MELCP1, MELCP2
  536. SEGDES MCHAM1
  537. C
  538. C
  539. MCHEL1=ICVN
  540. MCHAM1=MCHEL1.ICHAML(ISOUS)
  541. SEGACT MCHAM1
  542. MELCV1=MCHAM1.IELVAL(1)
  543. MELCV2=MCHAM1.IELVAL(2)
  544. SEGACT MELCV1, MELCV2
  545. SEGDES MCHAM1
  546. C
  547. C******* MELEME : maillage elementaire de ELTFA
  548. C
  549. MELEME=MLESPG.LECT(ISOUS)
  550. SEGACT MELEME
  551. NBEL1=MELEME.NUM(/2)
  552. NP1=MELEME.NUM(/1)
  553. C
  554. C******* IPT1 : SPG des inconnues primales et duales
  555. C
  556. IPT1=MLEPRI.LECT(ISOUS)
  557. NBEL = IPT1.NUM(/2)
  558. NP = IPT1.NUM(/1)
  559. C
  560. C******* MELGRR: SPG elementaire des coeffs de GRAD
  561. C Il resemble a IPT1, mais IPT1 peut changer!
  562. C
  563. MELGRR=MLEGRR.LECT(ISOUS)
  564. MELGRP=MLEGRP.LECT(ISOUS)
  565. MELGRV=MLEGRV.LECT(ISOUS)
  566. C
  567. C******* MLESPG et MLEPRI doivent avoir le meme nombre
  568. C d'elts et de noeuds!
  569. C
  570. IF((NBEL1 .NE. NBEL) .OR. (NP1 .NE. (NP-1)))THEN
  571. WRITE(IOIMP,*) 'ELTFA, SPGs gradients'
  572. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  573. CALL ERREUR(5)
  574. GOTO 9999
  575. ENDIF
  576. C
  577. MP = NP
  578. SEGINI RR , RUX , RUY , RRET ,
  579. & UXR , UXUX , UXUY , UXRET ,
  580. & UYR , UYUX , UYUY , UYRET ,
  581. & RETR , RETUX , RETUY , RETRET
  582. C
  583. C**** Duale = IMATRI.LISDUA(1) = 'RN'
  584. C Primale = IMATRI.LISPRI(1) = 'RN'
  585. C -> IMATRI.LIZAFM(1,1) = RR
  586. C
  587. C Duale = IMATRI.LISDUA(2) = 'RN'
  588. C Primale = IMATRI.LISPRI(2) = 'RUXN'
  589. C -> IMATRI.LIZAFM(1,2) = RUX
  590. C ...
  591. C
  592. C NB Pour le moment on y stoke le jacobien du flux
  593. C par rapport aux variables primitives (R,UX,UY,P)
  594. C
  595. IMATRI.LIZAFM(ISOUS,1) = RR
  596. IMATRI.LIZAFM(ISOUS,2) = RUX
  597. IMATRI.LIZAFM(ISOUS,3) = RUY
  598. IMATRI.LIZAFM(ISOUS,4) = RRET
  599. IMATRI.LIZAFM(ISOUS,5) = UXR
  600. IMATRI.LIZAFM(ISOUS,6) = UXUX
  601. IMATRI.LIZAFM(ISOUS,7) = UXUY
  602. IMATRI.LIZAFM(ISOUS,8) = UXRET
  603. IMATRI.LIZAFM(ISOUS,9) = UYR
  604. IMATRI.LIZAFM(ISOUS,10) = UYUX
  605. IMATRI.LIZAFM(ISOUS,11) = UYUY
  606. IMATRI.LIZAFM(ISOUS,12) = UYRET
  607. IMATRI.LIZAFM(ISOUS,13) = RETR
  608. IMATRI.LIZAFM(ISOUS,14) = RETUX
  609. IMATRI.LIZAFM(ISOUS,15) = RETUY
  610. IMATRI.LIZAFM(ISOUS,16) = RETRET
  611. C
  612. C******* Boucle sur les elts de MLEPRI.LECT
  613. C
  614.  
  615. DO IELPRI=1,NBEL,1
  616. IELEM=IELEM+1
  617. NGC=IPT1.NUM(NP,IELPRI)
  618. NLC=MLENTC.LECT(NGC)
  619. IF(NLC.NE.IELEM)THEN
  620. WRITE(IOIMP,*) 'NLC != IELEM'
  621. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  622. CALL ERREUR(5)
  623. GOTO 9999
  624. ENDIF
  625. VOLG = MPVOLU.VPOCHA(NLC,1)
  626. C
  627. C************* Pour la contribution des gradients
  628. C 1) les limiteurs
  629. C
  630. ALR = MPGLRN.VPOCHA(NLC,1)
  631. ALP = MPGLPN.VPOCHA(NLC,1)
  632. ALVX = MPGLVN.VPOCHA(NLC,1)
  633. ALVY = MPGLVN.VPOCHA(NLC,2)
  634. ICOOR= (NGC - 1) * 3
  635. XCG = MCOORD.XCOOR(ICOOR + 1)
  636. YCG = MCOORD.XCOOR(ICOOR + 2)
  637. C
  638. C******** Boucle sur le noeuds de MLESPG.LECT
  639. C
  640. DO IFAC=1,(NP - 1),1
  641. NGCF = MELEME.NUM(IFAC,IELPRI)
  642. ICOOR= (NGCF - 1) * 3
  643. DXCF = MCOORD.XCOOR(ICOOR + 1) - XCG
  644. DYCF = MCOORD.XCOOR(ICOOR + 2) - YCG
  645. NLCF = MLENTF.LECT(NGCF)
  646. NGCF1 = MELEFE.NUM(2,NLCF)
  647. IF(NGCF .NE. NGCF1)THEN
  648. WRITE(IOIMP,*)
  649. $ 'FACE ET FACEL = ???'
  650. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  651. CALL ERREUR(5)
  652. GOTO 9999
  653. ENDIF
  654. NGCG = MELEFE.NUM(1,NLCF)
  655. NGCD = MELEFE.NUM(3,NLCF)
  656. SURF = MPOVSU.VPOCHA(NLCF,1)
  657. NLFL = MLELIM.LECT(NGCF)
  658. IF(NLFL .NE. 0)THEN
  659. C
  660. C************* The point belongs on BC -> No contribution to jacobian!
  661. C
  662. DFRO(1)=0.0D0
  663. DFRO(2)=0.0D0
  664. DFRO(3)=0.0D0
  665. DFRO(4)=0.0D0
  666. DFRUN(1)=0.0D0
  667. DFRUN(2)=0.0D0
  668. DFRUN(3)=0.0D0
  669. DFRUN(4)=0.0D0
  670. DFRUT(1)=0.0D0
  671. DFRUT(2)=0.0D0
  672. DFRUT(3)=0.0D0
  673. DFRUT(4)=0.0D0
  674. DFRET(1)=0.0D0
  675. DFRET(2)=0.0D0
  676. DFRET(3)=0.0D0
  677. DFRET(4)=0.0D0
  678. C
  679. C************* On est sur G, sur D, ou sur un mur?
  680. C
  681. ELSEIF(NGCG .EQ. NGCD)THEN
  682. C
  683. C***************** Murs
  684. C Aucune contribution à l'elt droite!
  685. C
  686. CNX = MELVNX.VELCHE(1,NLCF)
  687. CNY = MELVNY.VELCHE(1,NLCF)
  688. CTX = MELT1X.VELCHE(1,NLCF)
  689. CTY = MELT1Y.VELCHE(1,NLCF)
  690. ROG = MELRO.VELCHE(1,NLCF)
  691. UNG = MELVUN.VELCHE(1,NLCF)
  692. UTG = MELVUT.VELCHE(1,NLCF)
  693. UXG = UNG * CNX + UTG * CTX
  694. UYG = UNG * CNY + UTG * CTY
  695. PG = MELP.VELCHE(1,NLCF)
  696. GAMG = MELGAM.VELCHE(1,NLCF)
  697. C
  698. C**************** Jacobien par rapport aux variables primitives sur
  699. C la face
  700. C
  701. IF(IIMPL .EQ. 4)THEN
  702. C
  703. C******************* VLH
  704. C
  705. CALL VLHJ3(ROG,UXG,UYG,PG,GAMG,CNX,CNY,DFRUN)
  706. C
  707. ELSEIF(IIMPL .EQ. 5)THEN
  708. C
  709. C******************* AUSM+
  710. C
  711. ROD=ROG
  712. PD=PG
  713. UND=-1.0D0*UNG
  714. UTD=UTG
  715. UXD = UND * CNX + UTD * CTX
  716. UYD = UND * CNY + UTD * CTY
  717. CALL AUSMW(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
  718. & GAMG,CNX,CNY,DFRUN)
  719. ELSE
  720. CALL ERREUR(5)
  721. ENDIF
  722. DFRO(1)=0.0D0
  723. DFRO(2)=0.0D0
  724. DFRO(3)=0.0D0
  725. DFRO(4)=0.0D0
  726. DFRUT(1)=0.0D0
  727. DFRUT(2)=0.0D0
  728. DFRUT(3)=0.0D0
  729. DFRUT(4)=0.0D0
  730. DFRET(1)=0.0D0
  731. DFRET(2)=0.0D0
  732. DFRET(3)=0.0D0
  733. DFRET(4)=0.0D0
  734. C
  735. ELSEIF(NGC .EQ. NGCG)THEN
  736. C
  737. C************* On est sur G
  738. C
  739. CNX = MELVNX.VELCHE(1,NLCF)
  740. CNY = MELVNY.VELCHE(1,NLCF)
  741. CTX = MELT1X.VELCHE(1,NLCF)
  742. CTY = MELT1Y.VELCHE(1,NLCF)
  743. C
  744. ROG = MELRO.VELCHE(1,NLCF)
  745. UNG = MELVUN.VELCHE(1,NLCF)
  746. UTG = MELVUT.VELCHE(1,NLCF)
  747. UXG = UNG * CNX + UTG * CTX
  748. UYG = UNG * CNY + UTG * CTY
  749. PG = MELP.VELCHE(1,NLCF)
  750. GAMG = MELGAM.VELCHE(1,NLCF)
  751. ECIN=0.5D0*ROG*(UNG*UNG+UTG*UTG)
  752. RETG=PG/(GAMG-1.0D0)+ECIN
  753. C
  754. ROD = MELRO.VELCHE(3,NLCF)
  755. UND = MELVUN.VELCHE(3,NLCF)
  756. UTD = MELVUT.VELCHE(3,NLCF)
  757. UXD = UND * CNX + UTD * CTX
  758. UYD = UND * CNY + UTD * CTY
  759. PD = MELP.VELCHE(3,NLCF)
  760. C
  761. NLCD=MLENTC.LECT(NGCD)
  762. VOLD=MPVOLU.VPOCHA(NLCD,1)
  763. IF(IIMPL .EQ. 4)THEN
  764. C
  765. C******************* VLH
  766. C
  767. CALL VLHJ1(ROG,UXG,UYG,PG,RETG,GAMG,CNX,CNY,CTX,CTY
  768. $ ,DFRO,DFRUN,DFRUT,DFRET)
  769. ELSEIF(IIMPL .EQ. 5)THEN
  770. C
  771. C******************* AUSM+
  772. C
  773. CALL AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
  774. & GAMG,CNX,CNY,CTX,CTY,
  775. & DFRO,DFRUN,DFRUT,DFRET)
  776. ELSE
  777. CALL ERREUR(5)
  778. ENDIF
  779. ELSEIF(NGC .EQ. NGCD)THEN
  780. CNX = -1.0D0 * MELVNX.VELCHE(1,NLCF)
  781. CNY = -1.0D0 * MELVNY.VELCHE(1,NLCF)
  782. CTX = -1.0D0 * MELT1X.VELCHE(1,NLCF)
  783. CTY = -1.0D0 * MELT1Y.VELCHE(1,NLCF)
  784. C
  785. ROG = MELRO.VELCHE(3,NLCF)
  786. UNG = MELVUN.VELCHE(3,NLCF)
  787. UTG = MELVUT.VELCHE(3,NLCF)
  788. UXG = -1.0D0 * (UNG * CNX + UTG * CTX)
  789. UYG = -1.0D0 * (UNG * CNY + UTG * CTY)
  790. PG = MELP.VELCHE(3,NLCF)
  791. GAMG = MELGAM.VELCHE(3,NLCF)
  792. ECIN=0.5D0*ROG*(UNG*UNG+UTG*UTG)
  793. RETG=PG/(GAMG-1.0D0)+ECIN
  794. C
  795. ROD = MELRO.VELCHE(1,NLCF)
  796. UND = MELVUN.VELCHE(1,NLCF)
  797. UTD = MELVUT.VELCHE(1,NLCF)
  798. UXD = -1.0D0 * (UND * CNX + UTD * CTX)
  799. UYD = -1.0D0 * (UND * CNY + UTD * CTY)
  800. PD = MELP.VELCHE(1,NLCF)
  801. C
  802. NLCD=MLENTC.LECT(NGCG)
  803. VOLD=MPVOLU.VPOCHA(NLCD,1)
  804. IF(IIMPL .EQ. 4)THEN
  805. C
  806. C******************* VLH
  807. C
  808. CALL VLHJ1(ROG,UXG,UYG,PG,RETG,GAMG,CNX,CNY,CTX,CTY
  809. $ ,DFRO,DFRUN,DFRUT,DFRET)
  810. ELSEIF(IIMPL .EQ. 5)THEN
  811. C
  812. C******************* AUSM+
  813. C
  814. CALL AUSMF(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
  815. & GAMG,CNX,CNY,CTX,CTY,
  816. & DFRO,DFRUN,DFRUT,DFRET)
  817. ELSE
  818. CALL ERREUR(5)
  819. ENDIF
  820. ELSE
  821. WRITE(IOIMP,*) 'FACEL = ???'
  822. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  823. CALL ERREUR(5)
  824. GOTO 9999
  825. ENDIF
  826. C
  827. DO I1=1,4
  828. DFRUX(I1)=DFRUN(I1)*CNX+DFRUT(I1)*CTX
  829. DFRUY(I1)=DFRUN(I1)*CNY+DFRUT(I1)*CTY
  830. ENDDO
  831. C
  832. NGNO=MELGRR.NUM(IFAC,IELPRI)
  833. NGN1=MELGRP.NUM(IFAC,IELPRI)
  834. NGN2=MELGRV.NUM(IFAC,IELPRI)
  835. C
  836. IF((NGN1 .NE. NGNO) .OR. (NGN2 .NE. NGNO))THEN
  837. WRITE(IOIMP,*) 'SPG COEFF GRAD =???'
  838. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  839. CALL ERREUR(5)
  840. GOTO 9999
  841. ENDIF
  842. C
  843. NLNO=MLENTC.LECT(NGNO)
  844. IF(NLNO .EQ. 0)THEN
  845. C
  846. C**************** MURS
  847. C
  848. IPT1.NUM(IFAC,IELPRI)=NGC
  849.  
  850. C Primal variable in IPT1.NUM(IFAC,IELPRI)=IPT1.NUM(NP,IELPRI)
  851. C Dual variable in IPT1.NUM(NP,IELPRI)
  852. C
  853. FUNCEL = -1.0D0 * SURF / VOLG
  854. UXR.AM(IELPRI,IFAC,NP)=UXR.AM(IELPRI,IFAC,NP)
  855. & +DFRUX(1)*FUNCEL
  856. UYR.AM(IELPRI,IFAC,NP)=UYR.AM(IELPRI,IFAC,NP)
  857. & +DFRUY(1)*FUNCEL
  858. UXUX.AM(IELPRI,IFAC,NP)=UXUX.AM(IELPRI,IFAC,NP)
  859. & +DFRUX(2)*FUNCEL
  860. UYUX.AM(IELPRI,IFAC,NP)=UYUX.AM(IELPRI,IFAC,NP)
  861. & +DFRUY(2)*FUNCEL
  862. UXUY.AM(IELPRI,IFAC,NP)=UXUY.AM(IELPRI,IFAC,NP)
  863. & +DFRUX(3)*FUNCEL
  864. UYUY.AM(IELPRI,IFAC,NP)=UYUY.AM(IELPRI,IFAC,NP)
  865. & +DFRUY(3)*FUNCEL
  866. UXRET.AM(IELPRI,IFAC,NP)=UXRET.AM(IELPRI,IFAC,NP)
  867. & +DFRUX(4)*FUNCEL
  868. UYRET.AM(IELPRI,IFAC,NP)=UYRET.AM(IELPRI,IFAC,NP)
  869. & +DFRUY(4)*FUNCEL
  870.  
  871. ELSE
  872. C
  873. C Primal variable in IPT1.NUM(NP,IELPRI)
  874. C Dual variable in IPT1.NUM(NP,IELPRI)
  875. C
  876. FUNCEL=DFRO(1)*SURF
  877. RR.AM(IELPRI,NP,NP)=RR.AM(IELPRI,NP,NP)-
  878. & FUNCEL/VOLG
  879. C
  880. C Primal variable in NP
  881. C Dual variable in IFAC
  882. C
  883. RR.AM(IELPRI,NP,IFAC)=RR.AM(IELPRI,NP,IFAC)+FUNCEL
  884. $ /VOLD
  885. C
  886. FUNCEL=DFRO(2)*SURF
  887. RUX.AM(IELPRI,NP,NP)=RUX.AM(IELPRI,NP,NP) -
  888. & FUNCEL/VOLG
  889. RUX.AM(IELPRI,NP,IFAC)=RUX.AM(IELPRI,NP,IFAC)+ FUNCEL
  890. $ /VOLD
  891. C
  892. FUNCEL=DFRO(3)*SURF
  893. RUY.AM(IELPRI,NP,NP)=RUY.AM(IELPRI,NP,NP) -
  894. & FUNCEL/VOLG
  895. RUY.AM(IELPRI,NP,IFAC)=RUY.AM(IELPRI,NP,IFAC)+ FUNCEL
  896. $ /VOLD
  897. C
  898. FUNCEL=DFRO(4)*SURF
  899. RRET.AM(IELPRI,NP,NP)=RRET.AM(IELPRI,NP,NP) -
  900. & FUNCEL/VOLG
  901. RRET.AM(IELPRI,NP,IFAC)= RRET.AM(IELPRI,NP,IFAC)
  902. $ +FUNCEL/VOLD
  903. C
  904. C
  905. FUNCEL=DFRUX(1)*SURF
  906. UXR.AM(IELPRI,NP,NP)=UXR.AM(IELPRI,NP,NP) -
  907. & FUNCEL/VOLG
  908. UXR.AM(IELPRI,NP,IFAC)=UXR.AM(IELPRI,NP,IFAC)+FUNCEL
  909. $ /VOLD
  910. C
  911. FUNCEL=DFRUX(2)*SURF
  912. UXUX.AM(IELPRI,NP,NP)=UXUX.AM(IELPRI,NP,NP) -
  913. & FUNCEL/VOLG
  914. UXUX.AM(IELPRI,NP,IFAC)=UXUX.AM(IELPRI,NP,IFAC)+
  915. $ FUNCEL/VOLD
  916. C
  917. FUNCEL=DFRUX(3)*SURF
  918. UXUY.AM(IELPRI,NP,NP)=UXUY.AM(IELPRI,NP,NP) -
  919. & FUNCEL/VOLG
  920. UXUY.AM(IELPRI,NP,IFAC)=UXUY.AM(IELPRI,NP,IFAC)+
  921. $ FUNCEL/VOLD
  922. C
  923. FUNCEL=DFRUX(4)*SURF
  924. UXRET.AM(IELPRI,NP,NP)=UXRET.AM(IELPRI,NP,NP) -
  925. & FUNCEL/VOLG
  926. UXRET.AM(IELPRI,NP,IFAC)=UXRET.AM(IELPRI,NP,IFAC)+
  927. $ FUNCEL/VOLD
  928. C
  929. C
  930. FUNCEL=DFRUY(1)*SURF
  931. UYR.AM(IELPRI,NP,NP)=UYR.AM(IELPRI,NP,NP) -
  932. & FUNCEL/VOLG
  933. UYR.AM(IELPRI,NP,IFAC)=UYR.AM(IELPRI,NP,IFAC)+FUNCEL
  934. $ /VOLD
  935. C
  936. FUNCEL=DFRUY(2)*SURF
  937. UYUX.AM(IELPRI,NP,NP)=UYUX.AM(IELPRI,NP,NP) -
  938. & FUNCEL/VOLG
  939. UYUX.AM(IELPRI,NP,IFAC)=UYUX.AM(IELPRI,NP,IFAC)+
  940. $ FUNCEL/VOLD
  941. C
  942. FUNCEL=DFRUY(3)*SURF
  943. UYUY.AM(IELPRI,NP,NP)=UYUY.AM(IELPRI,NP,NP) -
  944. & FUNCEL/VOLG
  945. UYUY.AM(IELPRI,NP,IFAC)=UYUY.AM(IELPRI,NP,IFAC)+
  946. $ FUNCEL/VOLD
  947. C
  948. FUNCEL=DFRUY(4)*SURF
  949. UYRET.AM(IELPRI,NP,NP)=UYRET.AM(IELPRI,NP,NP) -
  950. & FUNCEL/VOLG
  951. UYRET.AM(IELPRI,NP,IFAC)=UYRET.AM(IELPRI,NP,IFAC)+
  952. $ FUNCEL/VOLD
  953. C
  954. C
  955. FUNCEL=DFRET(1)*SURF
  956. RETR.AM(IELPRI,NP,NP)=RETR.AM(IELPRI,NP,NP) -
  957. & FUNCEL/VOLG
  958. RETR.AM(IELPRI,NP,IFAC)=RETR.AM(IELPRI,NP,IFAC)+FUNCEL
  959. $ /VOLD
  960. C
  961. FUNCEL=DFRET(2)*SURF
  962. RETUX.AM(IELPRI,NP,NP)=RETUX.AM(IELPRI,NP,NP) -
  963. & FUNCEL/VOLG
  964. RETUX.AM(IELPRI,NP,IFAC)=RETUX.AM(IELPRI,NP,IFAC)+
  965. $ FUNCEL/VOLD
  966. C
  967. FUNCEL=DFRET(3)*SURF
  968. RETUY.AM(IELPRI,NP,NP)=RETUY.AM(IELPRI,NP,NP) -
  969. & FUNCEL/VOLG
  970. RETUY.AM(IELPRI,NP,IFAC)=RETUY.AM(IELPRI,NP,IFAC)+
  971. $ FUNCEL/VOLD
  972. C
  973. FUNCEL=DFRET(4)*SURF
  974. RETRET.AM(IELPRI,NP,NP)=RETRET.AM(IELPRI,NP,NP) -
  975. & FUNCEL/VOLG
  976. RETRET.AM(IELPRI,NP,IFAC)=RETRET.AM(IELPRI,NP,IFAC)+
  977. $ FUNCEL/VOLD
  978. ENDIF
  979. C
  980. C************* Gradients computed with EULESCAL or EULEVECT
  981. C Computation of gradients contributions
  982.  
  983. DO IFAC2 = 1, (NP-1), 1
  984. C
  985. C**************** On controlle que les melemes de ELTFA et
  986. C le SPG du grad ont le meme ordre
  987. C
  988. NGFGR=MELEME.NUM(IFAC2,IELPRI)
  989. NLFGR=MLENTF.LECT(NGFGR)
  990. NG=MELEFE.NUM(1,NLFGR)
  991. ND=MELEFE.NUM(3,NLFGR)
  992. C
  993. IF(NG .EQ. ND)THEN
  994. C
  995. C******************* On est sur un mur: EULESCAL
  996. C
  997. DO IFACGR = 1, (NP -1) , 1
  998. C
  999. C********************** Boucle sur les elts de MELGRR pour trouver la position
  1000. C correspondante
  1001. C
  1002. NGGRA=MELGRR.NUM(IFACGR,IELPRI)
  1003. IF(NGFGR .EQ. NGGRA) GOTO 10
  1004. ENDDO
  1005. WRITE(IOIMP,*) 'ELTFA et SPG grad= ???'
  1006. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  1007. CALL ERREUR(5)
  1008. GOTO 9999
  1009. 10 CONTINUE
  1010. C
  1011. C******************* Les coeff qui nous interesse sont donc dans la
  1012. C position IFACGR du gradient
  1013. C
  1014. CELLR = ALR * (DXCF * MELCR1.VELCHE(IFACGR,IELPRI)
  1015. $ +DYCF * MELCR2.VELCHE(IFACGR,IELPRI))
  1016. CELLP = ALP * (DXCF * MELCP1.VELCHE(IFACGR,IELPRI)
  1017. $ + DYCF * MELCP2.VELCHE(IFACGR,IELPRI))
  1018. C
  1019. C******************* La vitesse est un peux plus difficile
  1020. C (HP EULEVECT ou vitesse imposé)
  1021. C
  1022. C CViVj contribution de vj sur le grad de vi
  1023. C
  1024. NLVIMP = MLEVL.LECT(NGGRA)
  1025. IF(NLVIMP .NE.0)THEN
  1026. C
  1027. C*********************** NGGRA spg de conditions limites -> no
  1028. C contribution
  1029. C
  1030. CVXVX=0.0D0
  1031. CVXVY=0.0D0
  1032. CVYVX=0.0D0
  1033. CVYVY=0.0D0
  1034. ELSE
  1035. CNX1=MPNORM.VPOCHA(NLFGR,1)
  1036. CNY1=MPNORM.VPOCHA(NLFGR,2)
  1037. CTX1=-1.0D0*CNY1
  1038. CTY1=CNX1
  1039. FUNCEL=(DXCF * MELCV1.VELCHE(IFACGR,IELPRI)+DYCF
  1040. $ * MELCV2.VELCHE(IFACGR,IELPRI))
  1041. CVXVX=FUNCEL*(CTX1*CTX1-CNX1*CNX1)*ALVX
  1042. CVXVY=FUNCEL*(CTX1*CTY1-CNX1*CNY1)*ALVX
  1043. CVYVX=FUNCEL*(CTX1*CTY1-CNX1*CNY1)*ALVY
  1044. CVYVY=FUNCEL*(CTY1*CTY1-CNY1*CNY1)*ALVY
  1045. ENDIF
  1046. C
  1047. ELSE
  1048. DO IFACGR = 1, (NP -1) , 1
  1049. C
  1050. C********************** Boucle sur les elts de MELGRR pour trouver la position
  1051. C correspondante
  1052. C
  1053. NGGRA=MELGRR.NUM(IFACGR,IELPRI)
  1054. IF((NGGRA .EQ. NG) .OR. (NGGRA .EQ. ND)) GOTO 20
  1055. ENDDO
  1056. WRITE(IOIMP,*) 'ELTFA et SPG grad= ???'
  1057. WRITE(IOIMP,*) 'SUBROUTINE konja5.eso'
  1058. CALL ERREUR(5)
  1059. GOTO 9999
  1060. 20 CONTINUE
  1061. CELLR = ALR * (DXCF * MELCR1.VELCHE(IFACGR,IELPRI)
  1062. $ +DYCF * MELCR2.VELCHE(IFACGR,IELPRI))
  1063. CELLP = ALP * (DXCF * MELCP1.VELCHE(IFACGR,IELPRI)
  1064. $ + DYCF * MELCP2.VELCHE(IFACGR,IELPRI))
  1065. CVXVX = ALVX * (DXCF * MELCV1.VELCHE(IFACGR,IELPRI)
  1066. $ + DYCF * MELCV2.VELCHE(IFACGR,IELPRI))
  1067. CVXVY=0.0D0
  1068. CVYVY = ALVY * (DXCF * MELCV1.VELCHE(IFACGR,IELPRI)
  1069. $ + DYCF * MELCV2.VELCHE(IFACGR,IELPRI))
  1070. CVYVX = 0.0D0
  1071.  
  1072. ENDIF
  1073. C
  1074. C**************** Dans le cas d'un mur (NLN0=0) , IFACGR = primale qui donne
  1075. C contribution en IELPRI,IFACGR,IFAC OU (IELPRI,IFACGR,NP)
  1076. C
  1077. C Dans le cas non mur (NLN0=0) , IFACGR = primale qui donne
  1078. C contribution en IELPRI,IFACGR,NP et en IELPRI,IFAC,NP
  1079. C
  1080. IF(NLNO .EQ. 0)THEN
  1081. C
  1082. C**************** MURS
  1083. C
  1084. FUNCEL = -1.0D0 * SURF / VOLG * CELLR
  1085. UXR.AM(IELPRI,IFACGR,NP)=UXR.AM(IELPRI,IFACGR,NP)
  1086. & +DFRUX(1)*FUNCEL
  1087. UYR.AM(IELPRI,IFACGR,NP)=UYR.AM(IELPRI,IFACGR,NP)
  1088. & +DFRUY(1)*FUNCEL
  1089. C
  1090. FUNCEL = -1.0D0 * SURF / VOLG
  1091. UXUX.AM(IELPRI,IFACGR,NP)=UXUX.AM(IELPRI,IFACGR,NP)
  1092. & +FUNCEL * (DFRUX(2)*CVXVX + DFRUX(3) *CVYVX)
  1093. UYUX.AM(IELPRI,IFACGR,NP)=UYUX.AM(IELPRI,IFACGR,NP)
  1094. & +FUNCEL * (DFRUY(2)*CVXVX + DFRUY(3) *CVYVX)
  1095. C
  1096. FUNCEL = -1.0D0 * SURF / VOLG
  1097. UXUY.AM(IELPRI,IFACGR,NP)=UXUY.AM(IELPRI,IFACGR,NP)
  1098. & +FUNCEL * (DFRUX(3)*CVYVY + DFRUX(2) *CVXVY)
  1099. UYUY.AM(IELPRI,IFACGR,NP)=UYUY.AM(IELPRI,IFACGR,NP)
  1100. & +FUNCEL * (DFRUY(3)*CVYVY + DFRUY(2) *CVXVY)
  1101.  
  1102. C
  1103. FUNCEL = -1.0D0 * SURF / VOLG * CELLP
  1104. UXRET.AM(IELPRI,IFACGR,NP)=UXRET.AM(IELPRI,IFACGR
  1105. $ ,NP)+DFRUX(4)*FUNCEL
  1106. UYRET.AM(IELPRI,IFACGR,NP)=UYRET.AM(IELPRI,IFACGR
  1107. $ ,NP)+DFRUY(4)*FUNCEL
  1108. ELSE
  1109. C
  1110. C******************* dR/
  1111. C
  1112. FUNCEL=DFRO(1)*SURF*CELLR
  1113. RR.AM(IELPRI,IFACGR,NP)=RR.AM(IELPRI,IFACGR,NP)-
  1114. & FUNCEL/VOLG
  1115. RR.AM(IELPRI,IFACGR,IFAC)=RR.AM(IELPRI,IFACGR,IFAC)
  1116. $ +FUNCEL/VOLD
  1117. C
  1118. FUNCEL=DFRO(4)*SURF*CELLP
  1119. RRET.AM(IELPRI,IFACGR,NP)=RRET.AM(IELPRI,IFACGR,NP)
  1120. $ -FUNCEL/VOLG
  1121. RRET.AM(IELPRI,IFACGR,IFAC)=RRET.AM(IELPRI,IFACGR
  1122. $ ,IFAC)+FUNCEL/VOLD
  1123. C
  1124. FUNCEL=SURF * (DFRO(2)*CVXVX + DFRO(3) *CVYVX)
  1125. RUX.AM(IELPRI,IFACGR,NP)=RUX.AM(IELPRI,IFACGR,NP)
  1126. & - FUNCEL / VOLG
  1127. RUX.AM(IELPRI,IFACGR,IFAC)=RUX.AM(IELPRI,IFACGR
  1128. $ ,IFAC) + FUNCEL / VOLD
  1129. C
  1130. FUNCEL=SURF * (DFRO(3)*CVYVY + DFRO(2) *CVXVY)
  1131. RUY.AM(IELPRI,IFACGR,NP)=RUY.AM(IELPRI,IFACGR,NP)
  1132. & - FUNCEL / VOLG
  1133. RUY.AM(IELPRI,IFACGR,IFAC)=RUY.AM(IELPRI,IFACGR
  1134. $ ,IFAC) + FUNCEL / VOLD
  1135. C
  1136. C******************* dUX/
  1137. C
  1138. FUNCEL=DFRUX(1)*SURF*CELLR
  1139. UXR.AM(IELPRI,IFACGR,NP)=UXR.AM(IELPRI,IFACGR,NP) -
  1140. & FUNCEL/VOLG
  1141. UXR.AM(IELPRI,IFACGR,IFAC)= UXR.AM(IELPRI,IFACGR
  1142. $ ,IFAC)+FUNCEL/VOLD
  1143. C
  1144. FUNCEL=DFRUX(4)*SURF*CELLP
  1145. UXRET.AM(IELPRI,IFACGR,NP)=UXRET.AM(IELPRI,IFACGR
  1146. $ ,NP) -FUNCEL/VOLG
  1147. UXRET.AM(IELPRI,IFACGR,IFAC)= UXRET.AM(IELPRI
  1148. $ ,IFACGR,IFAC)+FUNCEL/VOLD
  1149. C
  1150. FUNCEL=SURF * (DFRUX(2)*CVXVX + DFRUX(3) *CVYVX)
  1151. UXUX.AM(IELPRI,IFACGR,NP)=UXUX.AM(IELPRI,IFACGR,NP)
  1152. & - FUNCEL / VOLG
  1153. UXUX.AM(IELPRI,IFACGR,IFAC)=UXUX.AM(IELPRI,IFACGR
  1154. $ ,IFAC) + FUNCEL / VOLD
  1155. C
  1156. FUNCEL=SURF * (DFRUX(3)*CVYVY + DFRUX(2) *CVXVY)
  1157. UXUY.AM(IELPRI,IFACGR,NP)=UXUY.AM(IELPRI,IFACGR,NP)
  1158. & - FUNCEL / VOLG
  1159. UXUY.AM(IELPRI,IFACGR,IFAC)=UXUY.AM(IELPRI,IFACGR
  1160. $ ,IFAC) + FUNCEL / VOLD
  1161. C
  1162. C******************* dUY/
  1163. C
  1164. FUNCEL=DFRUY(1)*SURF*CELLR
  1165. UYR.AM(IELPRI,IFACGR,NP)=UYR.AM(IELPRI,IFACGR,NP) -
  1166. & FUNCEL/VOLG
  1167. UYR.AM(IELPRI,IFACGR,IFAC)=UYR.AM(IELPRI,IFACGR
  1168. $ ,IFAC) + FUNCEL/VOLD
  1169. C
  1170. FUNCEL=DFRUY(4)*SURF*CELLP
  1171. UYRET.AM(IELPRI,IFACGR,NP)=UYRET.AM(IELPRI,IFACGR
  1172. $ ,NP) -FUNCEL/VOLG
  1173. UYRET.AM(IELPRI,IFACGR,IFAC)=UYRET.AM(IELPRI,IFACGR
  1174. $ ,IFAC) + FUNCEL/VOLD
  1175. C
  1176. FUNCEL=SURF * (DFRUY(2)*CVXVX + DFRUY(3) *CVYVX)
  1177. UYUX.AM(IELPRI,IFACGR,NP)=UYUX.AM(IELPRI,IFACGR,NP)
  1178. & - FUNCEL / VOLG
  1179. UYUX.AM(IELPRI,IFACGR,IFAC)=UYUX.AM(IELPRI,IFACGR
  1180. $ ,IFAC) + FUNCEL / VOLD
  1181. C
  1182. FUNCEL=SURF * (DFRUY(3)*CVYVY + DFRUY(2) *CVXVY)
  1183. UYUY.AM(IELPRI,IFACGR,NP)=UYUY.AM(IELPRI,IFACGR,NP)
  1184. & - FUNCEL / VOLG
  1185. UYUY.AM(IELPRI,IFACGR,IFAC)=UYUY.AM(IELPRI,IFACGR
  1186. $ ,IFAC) + FUNCEL / VOLD
  1187. C
  1188. C******************* dRET/
  1189. C
  1190. FUNCEL=DFRET(1)*SURF*CELLR
  1191. RETR.AM(IELPRI,IFACGR,NP)=RETR.AM(IELPRI,IFACGR,NP)
  1192. $ -FUNCEL/VOLG
  1193. RETR.AM(IELPRI,IFACGR,IFAC)=RETR.AM(IELPRI,IFACGR
  1194. $ ,IFAC)+FUNCEL/VOLD
  1195. C
  1196. FUNCEL=DFRET(4)*SURF*CELLP
  1197. RETRET.AM(IELPRI,IFACGR,NP)=RETRET.AM(IELPRI,IFACGR
  1198. $ ,NP)-FUNCEL/VOLG
  1199. RETRET.AM(IELPRI,IFACGR,IFAC)=RETRET.AM(IELPRI
  1200. $ ,IFACGR,IFAC)+FUNCEL/VOLD
  1201. C
  1202. FUNCEL=SURF * (DFRET(2)*CVXVX + DFRET(3) *CVYVX)
  1203. RETUX.AM(IELPRI,IFACGR,NP)=RETUX.AM(IELPRI,IFACGR
  1204. $ ,NP)- FUNCEL / VOLG
  1205. RETUX.AM(IELPRI,IFACGR,IFAC)=RETUX.AM(IELPRI,IFACGR
  1206. $ ,IFAC) + FUNCEL / VOLD
  1207. C
  1208. FUNCEL=SURF * (DFRET(3)*CVYVY + DFRET(2) *CVXVY)
  1209. RETUY.AM(IELPRI,IFACGR,NP)=RETUY.AM(IELPRI,IFACGR
  1210. $ ,NP)- FUNCEL / VOLG
  1211. RETUY.AM(IELPRI,IFACGR,IFAC)=RETUY.AM(IELPRI,IFACGR
  1212. $ ,IFAC) + FUNCEL / VOLD
  1213. C
  1214. ENDIF
  1215. ENDDO
  1216. ENDDO
  1217. ENDDO
  1218. C
  1219. C
  1220. C******* Transformation:
  1221. C jacobien par rapport aux variables primitives ->
  1222. C jacobien par rapport aux variables conservatives
  1223. C Optimisation du programme:
  1224. C on pourrais faire ça avant
  1225. C
  1226. NBEL=RR.AM(/1)
  1227. NP=RR.AM(/2)
  1228. MP=RR.AM(/3)
  1229. C
  1230. C******* Boucle sur les variables primales
  1231. C
  1232. DO IELETR=1,NBEL,1
  1233. DO IPRITR=1,NP,1
  1234. NGPRI=IPT1.NUM(IPRITR,IELETR)
  1235. NLPRI=MLENTC.LECT(NGPRI)
  1236. ROG=MPOVRN.VPOCHA(NLPRI,1)
  1237. PG=MPOVPN.VPOCHA(NLPRI,1)
  1238. GAMG=MPOGAM.VPOCHA(NLPRI,1)
  1239. UXG=MPOVVN.VPOCHA(NLPRI,1)
  1240. UYG=MPOVVN.VPOCHA(NLPRI,2)
  1241. GAMGM1=GAMG-1.0D0
  1242. CELCO1=UXG/ROG
  1243. CELCO2=UYG/ROG
  1244. CELCO3=0.5D0*GAMGM1*(UXG*UXG+UYG*UYG)
  1245. CELCO4=GAMGM1*UXG
  1246. CELCO5=GAMGM1*UYG
  1247. DO IDUATR=1,MP,1
  1248. DDRO=RR.AM(IELETR,IPRITR,IDUATR)
  1249. DDUX=RUX.AM(IELETR,IPRITR,IDUATR)
  1250. DDUY=RUY.AM(IELETR,IPRITR,IDUATR)
  1251. DDP=RRET.AM(IELETR,IPRITR,IDUATR)
  1252. RR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX)
  1253. & -(CELCO2*DDUY)+(CELCO3*DDP)
  1254. RUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG -
  1255. & CELCO4*DDP
  1256. RUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG -
  1257. & CELCO5*DDP
  1258. RRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP
  1259. C
  1260. DDRO=UXR.AM(IELETR,IPRITR,IDUATR)
  1261. DDUX=UXUX.AM(IELETR,IPRITR,IDUATR)
  1262. DDUY=UXUY.AM(IELETR,IPRITR,IDUATR)
  1263. DDP=UXRET.AM(IELETR,IPRITR,IDUATR)
  1264. UXR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX)
  1265. & -(CELCO2*DDUY)+(CELCO3*DDP)
  1266. UXUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG -
  1267. & CELCO4*DDP
  1268. UXUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG -
  1269. & CELCO5*DDP
  1270. UXRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP
  1271. C
  1272. DDRO=UYR.AM(IELETR,IPRITR,IDUATR)
  1273. DDUX=UYUX.AM(IELETR,IPRITR,IDUATR)
  1274. DDUY=UYUY.AM(IELETR,IPRITR,IDUATR)
  1275. DDP=UYRET.AM(IELETR,IPRITR,IDUATR)
  1276. UYR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX)
  1277. & -(CELCO2*DDUY)+(CELCO3*DDP)
  1278. UYUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG -
  1279. & CELCO4*DDP
  1280. UYUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG -
  1281. & CELCO5*DDP
  1282. UYRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP
  1283. C
  1284. DDRO=RETR.AM(IELETR,IPRITR,IDUATR)
  1285. DDUX=RETUX.AM(IELETR,IPRITR,IDUATR)
  1286. DDUY=RETUY.AM(IELETR,IPRITR,IDUATR)
  1287. DDP=RETRET.AM(IELETR,IPRITR,IDUATR)
  1288. RETR.AM(IELETR,IPRITR,IDUATR)=DDRO-(CELCO1*DDUX)
  1289. & -(CELCO2*DDUY)+(CELCO3*DDP)
  1290. RETUX.AM(IELETR,IPRITR,IDUATR)=DDUX/ROG -
  1291. & CELCO4*DDP
  1292. RETUY.AM(IELETR,IPRITR,IDUATR)=DDUY/ROG -
  1293. & CELCO5*DDP
  1294. RETRET.AM(IELETR,IPRITR,IDUATR)=GAMGM1*DDP
  1295. ENDDO
  1296. ENDDO
  1297. ENDDO
  1298. C
  1299. MELEME = MLESPG.LECT(ISOUS)
  1300. SEGDES MELEME
  1301. MELEME = MLEPRI.LECT(ISOUS)
  1302. SEGDES MELEME
  1303. MELEME = MLEGRR.LECT(ISOUS)
  1304. SEGDES MELEME
  1305. MELEME = MLEGRP.LECT(ISOUS)
  1306. SEGDES MELEME
  1307. MELEME = MLEGRV.LECT(ISOUS)
  1308. SEGDES MELEME
  1309. C
  1310. SEGDES RR , RUX , RUY , RRET ,
  1311. & UXR , UXUX , UXUY , UXRET ,
  1312. & UYR , UYUX , UYUY , UYRET ,
  1313. & RETR , RETUX , RETUY , RETRET
  1314. SEGDES MELCR1, MELCR2
  1315. SEGDES MELCP1, MELCP2
  1316. SEGDES MELCV1, MELCV2
  1317. ENDDO
  1318. C
  1319. C
  1320. C**** Les MELVALs
  1321. C
  1322. SEGDES MELRO
  1323. SEGDES MELP
  1324. SEGDES MELGAM
  1325. SEGDES MELVUN
  1326. SEGDES MELVUT
  1327. SEGDES MELVNX
  1328. SEGDES MELVNY
  1329. SEGDES MELT1X
  1330. SEGDES MELT1Y
  1331. C
  1332. C**** Les MPOVALs
  1333. C
  1334. C
  1335. SEGDES MPOVRN
  1336. SEGDES MPOVVN
  1337. SEGDES MPOVPN
  1338. SEGDES MPOGAM
  1339. SEGDES MPVGRN
  1340. SEGDES MPGLRN
  1341. SEGDES MPVGPN
  1342. SEGDES MPGLPN
  1343. SEGDES MPVGPN
  1344. SEGDES MPGLPN
  1345. IF(MPVLIM .NE. 0) SEGDES MPVLIM
  1346. C
  1347. C**** Les MELEMEs, les MLENTIs, les volumes, les surfaces
  1348. C
  1349. SEGDES MELEMC
  1350. SEGSUP MLENTC
  1351. SEGSUP MLENTF
  1352. SEGSUP MLEVL
  1353. SEGDES MELEFE
  1354. SEGDES MPVOLU
  1355. SEGDES MPNORM
  1356. SEGDES MPOVSU
  1357. SEGDES MELPRI
  1358. SEGSUP MLEPRI
  1359. SEGSUP MLEGRR
  1360. SEGSUP MLEGRP
  1361. SEGSUP MLEGRV
  1362.  
  1363. C
  1364. C**** Coeff GRAD
  1365. C
  1366. MCHEL1=ICRN
  1367. SEGDES MCHEL1
  1368. MCHEL1=ICPN
  1369. SEGDES MCHEL1
  1370. MCHEL1=ICVN
  1371. SEGDES MCHEL1
  1372. C
  1373. SEGSUP MLESPG
  1374. SEGDES MELTFA
  1375. C
  1376. C**** MATRIK, liste des inconnues
  1377. C
  1378. SEGDES MATRIK
  1379. SEGDES IMATRI
  1380. SEGDES MLMINC
  1381.  
  1382. C
  1383. C
  1384. SEGDES RR , RUX , RUY , RRET ,
  1385. & UXR , UXUX , UXUY , UXRET ,
  1386. & UYR , UYUX , UYUY , UYRET ,
  1387. & RETR , RETUX , RETUY , RETRET
  1388.  
  1389. SEGSUP MLENTC
  1390. SEGSUP MLENTF
  1391. SEGDES MLMINC
  1392. SEGDES MLELIM
  1393.  
  1394. 9999 CONTINUE
  1395. RETURN
  1396. END
  1397.  
  1398.  
  1399.  
  1400.  
  1401.  
  1402.  
  1403.  
  1404.  
  1405.  
  1406.  
  1407.  

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