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

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