Télécharger konjr1.eso

Retour à la liste

Numérotation des lignes :

  1. C KONJR1 SOURCE PV 16/11/17 22:00:10 9180
  2. SUBROUTINE KONJR1(ILINC,IRN,IUN,IPN,IGAMN,INORM,ICHPVO
  3. $ ,ICHPSU,IUINF,IUPRI,MELEMC,MELEFE,MELLIM,IMAT)
  4. C
  5. C************************************************************************
  6. C
  7. C PROJET : CASTEM 2000
  8. C
  9. C NOM : KONJR1
  10. C
  11. C DESCRIPTION : Voir KON12
  12. C Calcul du jacobien du résidu pour la méthode
  13. C RUSANOLM (CENTERED + diff. num)
  14. C
  15. C Cas deux dimensions, gaz "calorically perfect"
  16. C
  17. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  18. C
  19. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF
  20. C
  21. C************************************************************************
  22. C
  23. C
  24. C APPELES (Outils
  25. C CASTEM) : KRIPAD, LICHT, ERREUR
  26. C
  27. C APPELES (Calcul) : CENTJ0, CENTJ2, FRBM1
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES
  32. C
  33. C ILINC : liste des inconnues (pointeur d'un objet de type LISTMOTS)
  34. C
  35. C 1) Pointeurs des CHPOINT
  36. C
  37. C IRN : CHPOINT CENTRE contenant la masse volumique ;
  38. C
  39. C IUN : CHPOINT CENTRE contenant la vitesse ;
  40. C
  41. C IPN : CHPOINT CENTRE contenant la pression ;
  42. C
  43. C IGAMN : CHPOINT CENTRE contenant le gamma ;
  44. C
  45. C INORM : CHPOINT FACE contenant les normales aux faces ;
  46. C
  47. C ICHPVO : CHPOINT VOLUME contenant le volume
  48. C
  49. C ICHPSU : CHPOINT FACE contenant la surface des faces
  50. C
  51. C
  52. C 2) Pointeurs de MELEME de la table DOMAINE
  53. C
  54. C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
  55. C
  56. C MELEFE : MELEME 'FACEL' du connectivité Faces -> Elts
  57. C
  58. C MELLIM : MELEME sur lequel on ne veut pas calculer
  59. C la contribution à la matrice jacobienne
  60. C
  61. C 3) Cas Bas Mach
  62. C
  63. C IUINF : CHPOINT, one component, "cut-off velocity"
  64. C 0 if there is no Bas MACH
  65. C
  66. C IUPRI : CHPOINT, one component, second "cut-off velocity"
  67. C 0 if there is no Bas MACH
  68. C
  69. C
  70. C SORTIES
  71. C
  72. C IMAT : pointeur de la MATRIK du jacobien du residu
  73. C
  74. C************************************************************************
  75. C
  76. C HISTORIQUE (Anomalies et modifications éventuelles)
  77. C
  78. C HISTORIQUE :
  79. C
  80. C************************************************************************
  81. C
  82. C
  83. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  84. C GAMMA \in (1,3)
  85. C Si non il faut le faire!!!
  86. C
  87. C************************************************************************
  88. C
  89. C
  90. C**** Variables de COOPTIO
  91. C
  92. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  93. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  94. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  95. C & ,IECHO, IIMPI, IOSPI
  96. C & ,IDIM
  97. C & ,MCOORD
  98. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  99. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  100. C & ,NORINC,NORVAL,NORIND,NORVAD
  101. C & ,NUCROU, IPSAUV
  102. C & ,IFICLE,IPREFI,IREFOR,ISAFOR
  103. C
  104. IMPLICIT INTEGER(I-N)
  105. INTEGER ILINC, IRN,IUN,IPN,IGAMN,INORM,ICHPVO,ICHPSU, IUINF, IUPRI
  106. & , IMAT, IGEOMC, IGEOMF
  107. & , NFAC, NBSOUS, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID
  108. & , NKMT, NBME, NBEL, MP, NP
  109. & , IFAC, NGCF, NLCF, NGCG, NGCD, NLCG, NLCD, NLFL, I1
  110. REAL*8 ROG, PG, UXG, UYG, GAMG, VOLG
  111. & , ROD, PD, UXD, UYD, GAMD, VOLD
  112. & , SURF, CNX, CNY, CTX, CTY, FUNCEL
  113. & , DFRO(4), DFRET(4), DFRUN(4), DFRUT(4), VINF
  114. & , UNG, UTG, UND, UTD, ROF, PF, UNF, UTF, GAMF, AMAT(4,4)
  115. & , RLAMMA, BMAT(4,4)
  116. CHARACTER*8 TYPE
  117. C
  118. C**** LES INCLUDES
  119. C
  120.  
  121. -INC PPARAM
  122. -INC CCOPTIO
  123. -INC SMCHPOI
  124. -INC SMELEME
  125. -INC SMLMOTS
  126. -INC SMLENTI
  127. POINTEUR MPRN.MPOVAL, MPUN.MPOVAL, MPPN.MPOVAL, MPGAMN.MPOVAL,
  128. & MPNORM.MPOVAL, MPVOLU.MPOVAL, MPOVSU.MPOVAL,
  129. & MPUPRI.MPOVAL, MPUINF.MPOVAL
  130. POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME,
  131. & MELEDU.MELEME, MELLIM.MELEME
  132. POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLELIM.MLENTI
  133. POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RRET.IZAFM,
  134. & UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXRET.IZAFM,
  135. & UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYRET.IZAFM,
  136. & RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETRET.IZAFM
  137. POINTEUR MLMINC.MLMOTS
  138. C
  139. C**** KRIPAD pour la correspondance global/local des conditions limits
  140. C
  141. CALL KRIPAD(MELLIM,MLELIM)
  142. c SEGACT MELLIM
  143. C
  144. C**** KRIPAD pour la correspondance global/local des centres
  145. C
  146. CALL KRIPAD(MELEMC,MLENTC)
  147. C
  148. C SEGACT MLENTC
  149. SEGACT MELEMC
  150. C
  151. SEGACT MELEFE
  152. C
  153. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  154. CALL LICHT(INORM,MPNORM,TYPE,IGEOMF)
  155. CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC)
  156. C
  157. C**** LICHT active les MPOVALs en *MOD
  158. C
  159. C i.e.
  160. C
  161. C SEGACT MPOVSU*MOD
  162. C SEGACT MPOVNO*MOD
  163. C SEGACT MPVOLU*MOD
  164. C
  165. MELEMF = IGEOMF
  166. CALL KRIPAD(MELEMF,MLENTF)
  167. C
  168. C SEGACT MLENTF
  169. SEGACT MELEMF
  170. C
  171. CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
  172. CALL LICHT(IPN,MPPN,TYPE,IGEOMC)
  173. CALL LICHT(IUN,MPUN,TYPE,IGEOMC)
  174. CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
  175. C
  176. C SEGACT MPRN*MOD
  177. C SEGACT MPPN*MOD
  178. C SEGACT MPUN*MOD
  179. C SEGACT MPGAMN*MOD
  180. C
  181. C**************************************************************
  182. C Bas Mach
  183. C**************************************************************
  184. CALL LICHT(IUPRI,MPUPRI,TYPE,IGEOMC)
  185. CALL LICHT(IUINF,MPUINF,TYPE,IGEOMC)
  186. C
  187. NFAC = MELEFE.NUM(/2)
  188. C
  189. C**** Maillage des inconnues primales
  190. C
  191. NBSOUS = 0
  192. NBREF = 0
  193. NBELEM = NFAC
  194. NBNN = 2
  195. C
  196. SEGINI MELEDU
  197. C MELEPR = MELEDU
  198. C
  199. C**** MELEDU = 'SEG2'
  200. C
  201. MELEDU.ITYPEL = 2
  202. C
  203. NRIGE = 7
  204. NMATRI = 1
  205. NKID = 9
  206. NKMT = 7
  207. C
  208. SEGINI MATRIK
  209. IMAT = MATRIK
  210. MATRIK.IRIGEL(1,1) = MELEDU
  211. MATRIK.IRIGEL(2,1) = MELEDU
  212. C
  213. C**** Matrice non symetrique
  214. C
  215. MATRIK.IRIGEL(7,1) = 2
  216. C
  217. NBME = 16
  218. NBSOUS = 1
  219. SEGINI IMATRI
  220. MLMINC = ILINC
  221. SEGACT MLMINC
  222. MATRIK.IRIGEL(4,1) = IMATRI
  223. C
  224. IMATRI.LISPRI(1) = MLMINC.MOTS(1)
  225. IMATRI.LISPRI(2) = MLMINC.MOTS(2)
  226. IMATRI.LISPRI(3) = MLMINC.MOTS(3)
  227. IMATRI.LISPRI(4) = MLMINC.MOTS(4)
  228. IMATRI.LISPRI(5) = MLMINC.MOTS(1)
  229. IMATRI.LISPRI(6) = MLMINC.MOTS(2)
  230. IMATRI.LISPRI(7) = MLMINC.MOTS(3)
  231. IMATRI.LISPRI(8) = MLMINC.MOTS(4)
  232. IMATRI.LISPRI(9) = MLMINC.MOTS(1)
  233. IMATRI.LISPRI(10) = MLMINC.MOTS(2)
  234. IMATRI.LISPRI(11) = MLMINC.MOTS(3)
  235. IMATRI.LISPRI(12) = MLMINC.MOTS(4)
  236. IMATRI.LISPRI(13) = MLMINC.MOTS(1)
  237. IMATRI.LISPRI(14) = MLMINC.MOTS(2)
  238. IMATRI.LISPRI(15) = MLMINC.MOTS(3)
  239. IMATRI.LISPRI(16) = MLMINC.MOTS(4)
  240. C
  241. IMATRI.LISDUA(1) = MLMINC.MOTS(1)
  242. IMATRI.LISDUA(2) = MLMINC.MOTS(1)
  243. IMATRI.LISDUA(3) = MLMINC.MOTS(1)
  244. IMATRI.LISDUA(4) = MLMINC.MOTS(1)
  245. IMATRI.LISDUA(5) = MLMINC.MOTS(2)
  246. IMATRI.LISDUA(6) = MLMINC.MOTS(2)
  247. IMATRI.LISDUA(7) = MLMINC.MOTS(2)
  248. IMATRI.LISDUA(8) = MLMINC.MOTS(2)
  249. IMATRI.LISDUA(9) = MLMINC.MOTS(3)
  250. IMATRI.LISDUA(10) = MLMINC.MOTS(3)
  251. IMATRI.LISDUA(11) = MLMINC.MOTS(3)
  252. IMATRI.LISDUA(12) = MLMINC.MOTS(3)
  253. IMATRI.LISDUA(13) = MLMINC.MOTS(4)
  254. IMATRI.LISDUA(14) = MLMINC.MOTS(4)
  255. IMATRI.LISDUA(15) = MLMINC.MOTS(4)
  256. IMATRI.LISDUA(16) = MLMINC.MOTS(4)
  257. C
  258. NBEL = NBELEM
  259. NBSOUS = 1
  260. NP = 2
  261. MP = 2
  262. SEGINI RR , RUX , RUY , RRET ,
  263. & UXR , UXUX , UXUY , UXRET ,
  264. & UYR , UYUX , UYUY , UYRET ,
  265. & RETR , RETUX , RETUY , RETRET
  266. C
  267. C**** Duale = IMATRI.LISDUA(1) = 'RN'
  268. C Primale = IMATRI.LISPRI(1) = 'RN'
  269. C -> IMATRI.LIZAFM(1,1) = RR
  270. C
  271. C Duale = IMATRI.LISDUA(2) = 'RN'
  272. C Primale = IMATRI.LISPRI(1) = 'RUXN'
  273. C -> IMATRI.LIZAFM(1,2) = RUX
  274. C ...
  275. C
  276. IMATRI.LIZAFM(1,1) = RR
  277. IMATRI.LIZAFM(1,2) = RUX
  278. IMATRI.LIZAFM(1,3) = RUY
  279. IMATRI.LIZAFM(1,4) = RRET
  280. IMATRI.LIZAFM(1,5) = UXR
  281. IMATRI.LIZAFM(1,6) = UXUX
  282. IMATRI.LIZAFM(1,7) = UXUY
  283. IMATRI.LIZAFM(1,8) = UXRET
  284. IMATRI.LIZAFM(1,9) = UYR
  285. IMATRI.LIZAFM(1,10) = UYUX
  286. IMATRI.LIZAFM(1,11) = UYUY
  287. IMATRI.LIZAFM(1,12) = UYRET
  288. IMATRI.LIZAFM(1,13) = RETR
  289. IMATRI.LIZAFM(1,14) = RETUX
  290. IMATRI.LIZAFM(1,15) = RETUY
  291. IMATRI.LIZAFM(1,16) = RETRET
  292. C
  293. DO IFAC = 1, NFAC, 1
  294. NGCF = MELEFE.NUM(2,IFAC)
  295. NLCF = MLENTF.LECT(NGCF)
  296. IF(NLCF .NE. IFAC)THEN
  297. WRITE(IOIMP,*) 'Il ne faut pas jouer avec la table domaine'
  298. CALL ERREUR(5)
  299. GOTO 9999
  300. ENDIF
  301. NGCG = MELEFE.NUM(1,IFAC)
  302. NGCD = MELEFE.NUM(3,IFAC)
  303. NLFL = MLELIM.LECT(NGCF)
  304. IF(NLFL .NE. 0)THEN
  305. C
  306. C********** The point belongs on BC -> No contribution to jacobian!
  307. C
  308. MELEDU.NUM(1,IFAC) = NGCG
  309. MELEDU.NUM(2,IFAC) = NGCD
  310.  
  311. ELSEIF(NGCG .NE. NGCD)THEN
  312. C
  313. C********** Les MELEMEs
  314. C
  315. MELEDU.NUM(1,IFAC) = NGCG
  316. MELEDU.NUM(2,IFAC) = NGCD
  317. C
  318. C********** Les etats G et D
  319. C
  320. NLCG = MLENTC.LECT(NGCG)
  321. NLCD = MLENTC.LECT(NGCD)
  322. C
  323. ROG = MPRN.VPOCHA(NLCG,1)
  324. PG = MPPN.VPOCHA(NLCG,1)
  325. UXG = MPUN.VPOCHA(NLCG,1)
  326. UYG = MPUN.VPOCHA(NLCG,2)
  327. GAMG = MPGAMN.VPOCHA(NLCG,1)
  328. VOLG = MPVOLU.VPOCHA(NLCG,1)
  329. C
  330. ROD = MPRN.VPOCHA(NLCD,1)
  331. PD = MPPN.VPOCHA(NLCD,1)
  332. UXD = MPUN.VPOCHA(NLCD,1)
  333. UYD = MPUN.VPOCHA(NLCD,2)
  334. GAMD = MPGAMN.VPOCHA(NLCD,1)
  335. VOLD = MPVOLU.VPOCHA(NLCD,1)
  336. C
  337. C********** La normale G->D
  338. C La tangente
  339. C
  340. SURF = MPOVSU.VPOCHA(NLCF,1)
  341. CNX = MPNORM.VPOCHA(NLCF,1)
  342. CNY = MPNORM.VPOCHA(NLCF,2)
  343. CTX = -1.0D0 * CNY
  344. CTY = CNX
  345. C
  346. C********** Cut off
  347. C
  348. VINF=MPUPRI.VPOCHA(NLCG,1)
  349. VINF=MAX(VINF,MPUPRI.VPOCHA(NLCD,1))
  350. VINF=MAX(VINF,MPUINF.VPOCHA(NLCG,1))
  351. VINF=MAX(VINF,MPUINF.VPOCHA(NLCD,1))
  352. C
  353. C********** L'etat moyen au centre
  354. C
  355. UNG=UXG*CNX+UYG*CNY
  356. UTG=UXG*CTX+UYG*CTY
  357. UND=UXD*CNX+UYD*CNY
  358. UTD=UXD*CTX+UYD*CTY
  359. UNF=0.5D0*(UNG+UND)
  360. UTF=0.5D0*(UTG+UTD)
  361. ROF=0.5D0*(ROG+ROD)
  362. PF=0.5D0*(PG+PD)
  363. GAMF=0.5D0*(GAMG+GAMD)
  364. C
  365. C********** We include in the cut-off UNG,UND,UTG,UTD in order to
  366. C avoid low diffusivity in stagnation regions
  367. C
  368. VINF=MAX(VINF,((UNG*UNG+UTG*UTG)**0.5D0))
  369. VINF=MAX(VINF,((UND*UND+UTD*UTD)**0.5D0))
  370. C
  371. C********** Calcule de la matrice de diffusivité pour le variables
  372. C conservatives.
  373. C Dans le repaire (n,t) local
  374. C F(UL,Ur) = 0.5 (F(UL)+F(UR)+AMAT*RLAMMA*(UL-UR))
  375. C
  376. CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA)
  377. C
  378. C Ici U=(\rho,\rho un, \rho ut, \rho et)
  379. C Il faut calculer la contribution sur F dans le repaire
  380. C (n,t) par rapport aux variables
  381. C U=(\rho,\rho ux, \rho uy, \rho et)
  382. C
  383. DO I1 = 1,4,1
  384. BMAT(I1,1) = AMAT(I1,1)
  385. BMAT(I1,2) = AMAT(I1,2) * CNX + AMAT(I1,3) * CTX
  386. BMAT(I1,3) = AMAT(I1,2) * CNY + AMAT(I1,3) * CTY
  387. BMAT(I1,4) = AMAT(I1,4)
  388. ENDDO
  389. C
  390. C********** La contribution de Gauche
  391. C
  392. C Partie centrée
  393. C
  394. CALL CENTJ0(ROG,UXG,UYG,PG,GAMG,CNX,CNY,CTX,CTY,
  395. & DFRO,DFRUN,DFRUT,DFRET)
  396. C
  397. C********** Centrée + Decentrée
  398. C
  399. C DFRO = derivative of the numerical approx. of (\rho un)
  400. C in the frame (x,y) (centered part)
  401. C BMAT(1,*) = derivative of the numerical approx. of (\rho un)
  402. C in the frame (x,y) (diff. part)
  403. C
  404. DO I1=1,4,1
  405. DFRO(I1) = DFRO(I1) + 0.5D0 * RLAMMA * BMAT(1,I1)
  406. DFRUN(I1) = DFRUN(I1) + 0.5D0 * RLAMMA * BMAT(2,I1)
  407. DFRUT(I1) = DFRUT(I1) + 0.5D0 * RLAMMA * BMAT(3,I1)
  408. DFRET(I1) = DFRET(I1) + 0.5D0 * RLAMMA * BMAT(4,I1)
  409. ENDDO
  410. C
  411. C
  412. C********** AB.AM(IFAC,IPRIM,IDUAL)
  413. C A = nom de l'inconnu duale (Ro,rUX,rUY,RET)
  414. C B = nom de l'inconnu primale (Ro,rUX,rUY,RET)
  415. C IPRIM = 1, 2 -> G, D
  416. C IDUAL = 1, 2 -> G, D
  417. C i.e.
  418. C A_IDUAL = AB.AM(IFAC,IPRIM,IDUAL) * B_IPRIM + ...
  419. C
  420. C
  421. C********** Dual RN
  422. C
  423. FUNCEL = SURF * DFRO(1)
  424. RR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  425. RR.AM(IFAC,1,2) = FUNCEL / VOLD
  426. C
  427. FUNCEL = SURF * DFRO(2)
  428. RUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  429. RUX.AM(IFAC,1,2) = FUNCEL / VOLD
  430. C
  431. FUNCEL = SURF * DFRO(3)
  432. RUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  433. RUY.AM(IFAC,1,2) = FUNCEL / VOLD
  434. C
  435. FUNCEL = SURF * DFRO(4)
  436. RRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  437. RRET.AM(IFAC,1,2) = FUNCEL / VOLD
  438. C
  439. C********** Dual RUXN
  440. C
  441. FUNCEL = SURF * (DFRUN(1) * CNX + DFRUT(1) * CTX)
  442. UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  443. UXR.AM(IFAC,1,2) = FUNCEL / VOLD
  444. C
  445. FUNCEL = SURF * (DFRUN(2) * CNX + DFRUT(2) * CTX)
  446. UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  447. UXUX.AM(IFAC,1,2) = FUNCEL / VOLD
  448. C
  449. FUNCEL = SURF * (DFRUN(3) * CNX + DFRUT(3) * CTX)
  450. UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  451. UXUY.AM(IFAC,1,2) = FUNCEL / VOLD
  452. C
  453. FUNCEL = SURF * (DFRUN(4) * CNX + DFRUT(4) * CTX)
  454. UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  455. UXRET.AM(IFAC,1,2) = FUNCEL / VOLD
  456. C
  457. C********** Dual RUYN
  458. C
  459. FUNCEL = SURF * (DFRUN(1) * CNY + DFRUT(1) * CTY)
  460. UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  461. UYR.AM(IFAC,1,2) = FUNCEL / VOLD
  462. C
  463. FUNCEL = SURF * (DFRUN(2) * CNY + DFRUT(2) * CTY)
  464. UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  465. UYUX.AM(IFAC,1,2) = FUNCEL / VOLD
  466. C
  467. FUNCEL = SURF * (DFRUN(3) * CNY + DFRUT(3) * CTY)
  468. UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  469. UYUY.AM(IFAC,1,2) = FUNCEL / VOLD
  470. C
  471. FUNCEL = SURF * (DFRUN(4) * CNY + DFRUT(4) * CTY)
  472. UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  473. UYRET.AM(IFAC,1,2) = FUNCEL / VOLD
  474. C
  475. C********** Dual RETN
  476. C
  477. FUNCEL = SURF * DFRET(1)
  478. RETR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  479. RETR.AM(IFAC,1,2) = FUNCEL / VOLD
  480. C
  481. FUNCEL = SURF * DFRET(2)
  482. RETUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  483. RETUX.AM(IFAC,1,2) = FUNCEL / VOLD
  484. C
  485. FUNCEL = SURF * DFRET(3)
  486. RETUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  487. RETUY.AM(IFAC,1,2) = FUNCEL / VOLD
  488. C
  489. FUNCEL = SURF * DFRET(4)
  490. RETRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  491. RETRET.AM(IFAC,1,2) = FUNCEL / VOLD
  492. C
  493. C
  494. C********** La contribution de D
  495. C Je change l'orientation de normales
  496. C Dans ce cas F(UD,UG) = 0.5 * (F(UD)+F(UD)
  497. C + RLAMMA * AMAT * (UD - UG))
  498. C
  499. CNX = -1.0D0 * CNX
  500. CNY = -1.0D0 * CNY
  501. CTX = -1.0D0 * CTX
  502. CTY = -1.0D0 * CTY
  503. C
  504. C********** L'etat moyen au centre
  505. C Meme que avant; mais UNF, UTF change de signe
  506. C
  507. UNF=-1.0D0*UNF
  508. UTF=-1.0D0*UTF
  509. C
  510. C********** On recalcule de la matrice de diffusivité pour le variables
  511. C conservatives (pas d'optimisation illisible!!!)
  512. C
  513. CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA)
  514. C
  515. C Ici U=(\rho,\rho un, \rho ut, \rho et)
  516. C Il faut calculer la contribution sur F dans le repaire
  517. C (n,t) par rapport aux variables
  518. C U=(\rho,\rho ux, \rho uy, \rho et)
  519. C
  520. DO I1 = 1,4,1
  521. BMAT(I1,1) = AMAT(I1,1)
  522. BMAT(I1,2) = AMAT(I1,2) * CNX + AMAT(I1,3) * CTX
  523. BMAT(I1,3) = AMAT(I1,2) * CNY + AMAT(I1,3) * CTY
  524. BMAT(I1,4) = AMAT(I1,4)
  525. ENDDO
  526. C
  527. CALL CENTJ0(ROD,UXD,UYD,PD,GAMD,CNX,CNY,CTX,CTY,
  528. & DFRO,DFRUN,DFRUT,DFRET)
  529. C
  530. C********** Centrée + Decentrée
  531. C
  532. DO I1=1,4,1
  533. DFRO(I1) = DFRO(I1) + 0.5D0 * RLAMMA * BMAT(1,I1)
  534. DFRUN(I1) = DFRUN(I1) + 0.5D0 * RLAMMA * BMAT(2,I1)
  535. DFRUT(I1) = DFRUT(I1) + 0.5D0 * RLAMMA * BMAT(3,I1)
  536. DFRET(I1) = DFRET(I1) + 0.5D0 * RLAMMA * BMAT(4,I1)
  537. ENDDO
  538. C
  539. C
  540. C********** Dual RN
  541. C
  542. FUNCEL = SURF * DFRO(1)
  543. RR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  544. RR.AM(IFAC,2,1) = FUNCEL / VOLG
  545. C
  546. FUNCEL = SURF * DFRO(2)
  547. RUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  548. RUX.AM(IFAC,2,1) = FUNCEL / VOLG
  549. C
  550. FUNCEL = SURF * DFRO(3)
  551. RUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  552. RUY.AM(IFAC,2,1) = FUNCEL / VOLG
  553. C
  554. FUNCEL = SURF * DFRO(4)
  555. RRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  556. RRET.AM(IFAC,2,1) = FUNCEL / VOLG
  557. C
  558. C********** Dual RUXN
  559. C
  560. FUNCEL = SURF * (DFRUN(1) * CNX + DFRUT(1) * CTX)
  561. UXR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  562. UXR.AM(IFAC,2,1) = FUNCEL / VOLG
  563. C
  564. FUNCEL = SURF * (DFRUN(2) * CNX + DFRUT(2) * CTX)
  565. UXUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  566. UXUX.AM(IFAC,2,1) = FUNCEL / VOLG
  567. C
  568. FUNCEL = SURF * (DFRUN(3) * CNX + DFRUT(3) * CTX)
  569. UXUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  570. UXUY.AM(IFAC,2,1) = FUNCEL / VOLG
  571. C
  572. FUNCEL = SURF * (DFRUN(4) * CNX + DFRUT(4) * CTX)
  573. UXRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  574. UXRET.AM(IFAC,2,1) = FUNCEL / VOLG
  575. C
  576. C********** Dual RUYN
  577. C
  578. FUNCEL = SURF * (DFRUN(1) * CNY + DFRUT(1) * CTY)
  579. UYR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  580. UYR.AM(IFAC,2,1) = FUNCEL / VOLG
  581. C
  582. FUNCEL = SURF * (DFRUN(2) * CNY + DFRUT(2) * CTY)
  583. UYUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  584. UYUX.AM(IFAC,2,1) = FUNCEL / VOLG
  585. C
  586. FUNCEL = SURF * (DFRUN(3) * CNY + DFRUT(3) * CTY)
  587. UYUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  588. UYUY.AM(IFAC,2,1) = FUNCEL / VOLG
  589. C
  590. FUNCEL = SURF * (DFRUN(4) * CNY + DFRUT(4) * CTY)
  591. UYRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  592. UYRET.AM(IFAC,2,1) = FUNCEL / VOLG
  593. C
  594. C********** Dual RETN
  595. C
  596. FUNCEL = SURF * DFRET(1)
  597. RETR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  598. RETR.AM(IFAC,2,1) = FUNCEL / VOLG
  599. C
  600. FUNCEL = SURF * DFRET(2)
  601. RETUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  602. RETUX.AM(IFAC,2,1) = FUNCEL / VOLG
  603. C
  604. FUNCEL = SURF * DFRET(3)
  605. RETUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  606. RETUY.AM(IFAC,2,1) = FUNCEL / VOLG
  607. C
  608. FUNCEL = SURF * DFRET(4)
  609. RETRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
  610. RETRET.AM(IFAC,2,1) = FUNCEL / VOLG
  611. C
  612. ELSE
  613. C
  614. C********** Murs (NGCG = NGCD)
  615. C
  616. C
  617. C********** Les MELEMEs
  618. C
  619. MELEDU.NUM(1,IFAC) = NGCG
  620. MELEDU.NUM(2,IFAC) = NGCD
  621. NLCG = MLENTC.LECT(NGCG)
  622. C
  623. ROG = MPRN.VPOCHA(NLCG,1)
  624. PG = MPPN.VPOCHA(NLCG,1)
  625. UXG = MPUN.VPOCHA(NLCG,1)
  626. UYG = MPUN.VPOCHA(NLCG,2)
  627. GAMG = MPGAMN.VPOCHA(NLCG,1)
  628. VOLG = MPVOLU.VPOCHA(NLCG,1)
  629. C
  630. C********** Cut off
  631. C
  632. VINF=MPUPRI.VPOCHA(NLCG,1)
  633. VINF=MAX(VINF,MPUINF.VPOCHA(NLCG,1))
  634. C
  635. C********** La normale sortante
  636. C
  637. SURF = MPOVSU.VPOCHA(NLCF,1)
  638. CNX = MPNORM.VPOCHA(NLCF,1)
  639. CNY = MPNORM.VPOCHA(NLCF,2)
  640. C
  641. C********** L'etat moyen au centre
  642. C
  643. UNF=0.0D0
  644. UTF=0.5D0*(UXG*CTX+UYG*CTY)
  645. ROF=ROG
  646. PF=PG
  647. GAMF=GAMG
  648. C
  649. C********** We include in the cut-off UG
  650. C
  651. VINF=MAX(VINF,((UXG*UXG+UYG*UYG)**0.5D0))
  652. C
  653. C********** Calcule de la matrice de diffusivité pour le variables
  654. C conservatives.
  655. C Dans le repaire (n,t) local
  656. C F(UL,Ur) = 0.5 (F(UL)+F(UR)+AMAT*RLAMMA*(UL-UR))
  657. C
  658. CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA)
  659. C
  660. C Ici U=(\rho,\rho un, \rho ut, \rho et)
  661. C Il faut calculer la contribution sur F dans le repaire
  662. C (n,t) par rapport aux variables
  663. C U=(\rho,\rho ux, \rho uy, \rho et)
  664. C Dans le cas murs, le seul flux normal donne une contribution
  665. C non nulle
  666. C En plus AMAT(1,2)=AMAT(3,2)=AMAT(4,2)=0
  667. C AMAT(2,2)=1.0D0
  668. C
  669. C write(*,*) 'AMAT(*,2)=',(AMAT(I1,2),I1=1,4)
  670. C
  671. C********** Contribution centrée
  672. C
  673. CALL CENTJ2(ROG,UXG,UYG,GAMG,CNX,CNY,
  674. & DFRUN)
  675. C
  676. C********** Centrée + Decentrée
  677. C
  678. DFRUN(2) = DFRUN(2) + (RLAMMA * CNX)
  679. DFRUN(3) = DFRUN(3) + (RLAMMA * CNY)
  680. C
  681. C********** Dual RN
  682. C
  683. RR.AM(IFAC,1,1) = 0.0D0
  684. RR.AM(IFAC,1,2) = 0.0D0
  685. C
  686. RUX.AM(IFAC,1,1) = 0.0D0
  687. RUX.AM(IFAC,1,2) = 0.0D0
  688. C
  689. RUY.AM(IFAC,1,1) = 0.0D0
  690. RUY.AM(IFAC,1,2) = 0.0D0
  691. C
  692. RRET.AM(IFAC,1,1) = 0.0D0
  693. RRET.AM(IFAC,1,2) = 0.0D0
  694. C
  695. C********** Dual RUXN
  696. C
  697. FUNCEL = SURF * DFRUN(1) * CNX
  698. UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  699. UXR.AM(IFAC,1,2) = 0.0D0
  700. C
  701. FUNCEL = SURF * DFRUN(2) * CNX
  702. UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  703. UXUX.AM(IFAC,1,2) = 0.0D0
  704. C
  705. FUNCEL = SURF * DFRUN(3) * CNX
  706. UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  707. UXUY.AM(IFAC,1,2) = 0.0D0
  708. C
  709. FUNCEL = SURF * DFRUN(4) * CNX
  710. UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  711. UXRET.AM(IFAC,1,2) = 0.0D0
  712. C
  713. C********** Dual RUYN
  714. C
  715. FUNCEL = SURF * DFRUN(1) * CNY
  716. UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  717. UYR.AM(IFAC,1,2) = 0.0D0
  718. C
  719. FUNCEL = SURF * DFRUN(2) * CNY
  720. UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  721. UYUX.AM(IFAC,1,2) = 0.0D0
  722. C
  723. FUNCEL = SURF * DFRUN(3) * CNY
  724. UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  725. UYUY.AM(IFAC,1,2) = 0.0D0
  726. C
  727. FUNCEL = SURF * DFRUN(4) * CNY
  728. UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
  729. UYRET.AM(IFAC,1,2) = 0.0D0
  730. C
  731. C********** Dual RETN
  732. C
  733. RETR.AM(IFAC,1,1) = 0.0D0
  734. RETR.AM(IFAC,1,2) = 0.0D0
  735. C
  736. RETUX.AM(IFAC,1,1) = 0.0D0
  737. RETUX.AM(IFAC,1,2) = 0.0D0
  738. C
  739. RETUY.AM(IFAC,1,1) = 0.0D0
  740. RETUY.AM(IFAC,1,2) = 0.0D0
  741. C
  742. RETRET.AM(IFAC,1,1) = 0.0D0
  743. RETRET.AM(IFAC,1,2) = 0.0D0
  744. C
  745. C********** Dual RN
  746. C
  747. RR.AM(IFAC,2,2) = 0.0D0
  748. RR.AM(IFAC,2,1) = 0.0D0
  749. C
  750. RUX.AM(IFAC,2,2) = 0.0D0
  751. RUX.AM(IFAC,2,1) = 0.0D0
  752. C
  753. RUY.AM(IFAC,2,2) = 0.0D0
  754. RUY.AM(IFAC,2,1) = 0.0D0
  755. C
  756. RRET.AM(IFAC,2,2) = 0.0D0
  757. RRET.AM(IFAC,2,1) = 0.0D0
  758. C
  759. C********** Dual RUXN
  760. C
  761. UXR.AM(IFAC,2,2) = 0.0D0
  762. UXR.AM(IFAC,2,1) = 0.0D0
  763. C
  764. UXUX.AM(IFAC,2,2) = 0.0D0
  765. UXUX.AM(IFAC,2,1) = 0.0D0
  766. C
  767. UXUY.AM(IFAC,2,2) = 0.0D0
  768. UXUY.AM(IFAC,2,1) = 0.0D0
  769. C
  770. UXRET.AM(IFAC,2,2) = 0.0D0
  771. UXRET.AM(IFAC,2,1) = 0.0D0
  772. C
  773. C********** Dual RUYN
  774. C
  775. UYR.AM(IFAC,2,2) = 0.0D0
  776. UYR.AM(IFAC,2,1) = 0.0D0
  777. C
  778. UYUX.AM(IFAC,2,2) = 0.0D0
  779. UYUX.AM(IFAC,2,1) = 0.0D0
  780. C
  781. UYUY.AM(IFAC,2,2) = 0.0D0
  782. UYUY.AM(IFAC,2,1) = 0.0D0
  783. C
  784. UYRET.AM(IFAC,2,2) = 0.0D0
  785. UYRET.AM(IFAC,2,1) = 0.0D0
  786. C
  787. C********** Dual RETN
  788. C
  789. RETR.AM(IFAC,2,2) = 0.0D0
  790. RETR.AM(IFAC,2,1) = 0.0D0
  791. C
  792. RETUX.AM(IFAC,2,2) = 0.0D0
  793. RETUX.AM(IFAC,2,1) = 0.0D0
  794. C
  795. RETUY.AM(IFAC,2,2) = 0.0D0
  796. RETUY.AM(IFAC,2,1) = 0.0D0
  797. C
  798. RETRET.AM(IFAC,2,2) = 0.0D0
  799. RETRET.AM(IFAC,2,1) = 0.0D0
  800. C
  801. ENDIF
  802. ENDDO
  803. C
  804. SEGDES MELEMC
  805. SEGDES MELEFE
  806. SEGDES MELEMF
  807. C
  808. SEGDES MPOVSU
  809. SEGDES MPVOLU
  810. SEGDES MPNORM
  811. C
  812. SEGDES MPRN
  813. SEGDES MPPN
  814. SEGDES MPUN
  815. SEGDES MPGAMN
  816. C
  817. SEGDES MELEDU
  818. SEGDES MATRIK
  819. SEGDES IMATRI
  820. C
  821. SEGDES RR , RUX , RUY , RRET ,
  822. & UXR , UXUX , UXUY , UXRET ,
  823. & UYR , UYUX , UYUY , UYRET ,
  824. & RETR , RETUX , RETUY , RETRET
  825.  
  826. SEGSUP MLENTC
  827. SEGSUP MLENTF
  828. SEGDES MLMINC
  829. SEGSUP MLELIM
  830. C
  831. IF(MELLIM .NE.0) SEGDES MELLIM
  832. C
  833. SEGDES MPUPRI
  834. SEGDES MPUINF
  835. C
  836. 9999 CONTINUE
  837. RETURN
  838. END
  839.  
  840.  
  841.  
  842.  
  843.  
  844.  
  845.  
  846.  

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