Télécharger kfm13.eso

Retour à la liste

Numérotation des lignes :

  1. C KFM13 SOURCE PV 20/03/30 21:20:32 10567
  2. C KFM13 SOURCE CHAT 05/01/13 00:55:24 5004
  3. SUBROUTINE KFM13(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
  4. & ICHPSU,ICHPDI,ICHPVO,INORM,
  5. & MELEMC,MELEMF,MELEFE,DTPS,DCFL,
  6. & MELLIM,ICHLIM,NJAC,ICOEV,
  7. & IDU,IPROBL)
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : KFM13
  13. C
  14. C DESCRIPTION : Voir aussi KFM1
  15. C
  16. C Cas 3D, gaz "calorically perfect"
  17. C Méthode implicite sans Matrice résolue
  18. C par Jacobi par point.
  19. C
  20. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  21. C
  22. C AUTEUR : T. KLOCZKO, DEN/DM2S/SFME/LTMF
  23. C
  24. C************************************************************************
  25. C ENTREES
  26. C
  27. C
  28. C 1) Pointeurs des CHPOINTs
  29. C
  30. C IRES : CHPOINT 'CENTRE' contenant le residu
  31. C
  32. C IRN : CHPOINT 'CENTRE' contenant la masse volumique
  33. C
  34. C IGN : CHPOINT 'CENTRE' contenant la q.d.m.
  35. C
  36. C IRETN : CHPOINT 'CENTRE' contenant l'energie totale
  37. C
  38. C IGAMN : CHPOINT 'CENTRE' contenant le gamma
  39. C
  40. C IVCO : CHPOINT 'CENTRE' contenant le cut-off
  41. C
  42. C ICHLIM : CHPOINT conditions aux bords
  43. C
  44. C ICOEV : CHPOINT 'CENTRE' contenant le (gamma - 1) K / (R \rho)
  45. C coeff pour calculer le ray. spect. de la matrice
  46. C visc.
  47. C
  48. C 3) Pointeurs de CHPOINTs de la table DOMAINE
  49. C
  50. C ICHPSU : CHPOINT "FACE" contenant la surface des faces
  51. C
  52. C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
  53. C de chaque element
  54. C
  55. C ICHPVO : CHPOINT "CENTRE" contenant le volume
  56. C de chaque element
  57. C
  58. C INORM : CHPOINT "CENTRE" contenant le normales aux faces
  59. C
  60. C
  61. C 4) Pointeurs de MELEME de la table DOMAINE
  62. C
  63. C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
  64. C
  65. C MELEMF : MELEME 'FACE' du SPG des FACES
  66. C
  67. C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
  68. C
  69. C MELLIM : MELEME SPG conditions aux bords
  70. C
  71. C 5)
  72. C
  73. C DCFL = le double de la CFL
  74. C
  75. C DTPS = pas de temps physique
  76. C
  77. C NJAC = nombre d'iterations dans le jacobien
  78. C
  79. C
  80. C SORTIES
  81. C
  82. C IDU : resultat
  83. C
  84. C IPROBL : si > 0, il y a un probleme
  85. C
  86. C************************************************************************
  87. C
  88. C HISTORIQUE (Anomalies et modifications éventuelles)
  89. C
  90. C HISTORIQUE : 2004 : création
  91. C
  92. C************************************************************************
  93. C
  94. C
  95. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  96. C GAMMA \in (1,3)
  97. C Y \in (0,1)
  98. C Si non il faut le faire!!!
  99. C
  100. C************************************************************************
  101. C
  102. C
  103. C**** Variables de COOPTIO
  104. C
  105. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  106. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  107. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  108. C & ,IECHO, IIMPI, IOSPI
  109. C & ,IDIM
  110. CC & ,MCOORD
  111. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  112. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  113. C & ,NORINC,NORVAL,NORIND,NORVAD
  114. C & ,NUCROU, IPSAUV, IFICLE, IPREFI, IREFOR, ISAFOR
  115. C
  116. IMPLICIT INTEGER(I-N)
  117. INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPDI,ICHPVO,INORM
  118. & ,MELEMC,MELEMF,MELEFE,IDU,NFAC
  119. & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED
  120. & ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP
  121. & ,I1,I2, IBLJAC, IPROBL, IVCO, ICOEV
  122. REAL*8 ROG, RUXG, RUYG, RUZG, RETG, GAMG, RHTG, REC
  123. & , ROD, RUXD, RUYD, RUZD, RETD, GAMD, RHTD
  124. & , UXG, UYG, UZG, UNG, UTG, UVG, PG, VOLG, DIAMG, SURF
  125. & , UXD, UYD, UZD, UND, UTD, UVD, PD, VOLD, DIAMD
  126. & , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ
  127. & , DCFL, DTPS, UNSDT, UNSDTF
  128. & , FR, FRUX, FRUY, FRUZ, FRET
  129. & , QG, CG, UCOG, ROF, PF, UNF, GAMF, UCOF, CF, SIGMAF
  130. & , QN, PHI2, COEF, COEF1
  131. CHARACTER*8 TYPE
  132. C
  133. C**** LES INCLUDES
  134. C
  135.  
  136. -INC PPARAM
  137. -INC CCOPTIO
  138. -INC SMCHPOI
  139. POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
  140. & , MPODU.MPOVAL, MPOVOL.MPOVAL
  141. & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL
  142. & , MPNORM.MPOVAL, MPLIM.MPOVAL, MPORES.MPOVAL, MPOCO.MPOVAL
  143. & , MPOCO1.MPOVAL, MPCOEV.MPOVAL
  144. -INC SMCOORD
  145. -INC SMELEME
  146. -INC SMLMOTS
  147. -INC SMLENTI
  148. POINTEUR MLELIM.MLENTI
  149. C
  150. C**** POINTEUR vecteurs
  151. C
  152. SEGMENT VEC
  153. REAL*8 VALVEC(I1)
  154. ENDSEGMENT
  155. POINTEUR D1.VEC,SIGMA0.VEC, AN0.VEC, DN0.VEC
  156. & ,USDTI.VEC, PHI2V.VEC, SIGMAV.VEC, NEWLI.VEC
  157. C
  158. C**** POINTEUR matriciels
  159. C
  160. SEGMENT MAT
  161. REAL*8 VAL(I1,I2)
  162. ENDSEGMENT
  163. POINTEUR CONS1.MAT, CONS2.MAT, BN0.MAT, BN1.MAT
  164. & ,CN0.MAT, CN1.MAT, DU.MAT
  165. & , RES.MAT
  166. & , Q1.MAT, Q2.MAT
  167. C
  168. C**** IPROBL: si different de 0, un probleme a été detecté
  169. C
  170. IPROBL=0
  171. C
  172. C**** Initialisation des MLENTI des conditions aux limites
  173. C
  174. C
  175. CALL KRIPAD(MELLIM,MLELIM)
  176. C SEGINI MLELIM
  177. C
  178. C**** Initialisation des MELEMEs
  179. C
  180. C 'CENTRE', 'FACEL'
  181. C
  182. IPT2 = MELEFE
  183. SEGACT IPT2
  184. NFAC = IPT2.NUM(/2)
  185. C
  186. C**** KRIPAD pour la correspondance global/local de centre
  187. C
  188. CALL KRIPAD(MELEMC,MLENT1)
  189. C
  190. C**** MLENTI1 a nbpts elements
  191. C
  192. C Si i est le numero global d'un noeud de ICEN,
  193. C MLENT1.LECT(i) contient sa position, i.e.
  194. C
  195. C I = numero global du noeud centre
  196. C MLENT1.LECT(i) = numero local du noeud centre
  197. C
  198. C MLENT1 déjà activé, i.e.
  199. C
  200. C SEGINI MLENT1
  201. C
  202. C
  203. C**** KRIPAD pour la correspondance global/local de 'FACE'
  204. C
  205. CALL KRIPAD(MELEMF,MLENT2)
  206. C SEGINI MLENT2
  207. C
  208. C
  209. C**** CHPOINTs de la table DOMAINE
  210. C
  211. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  212. CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  213. CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC)
  214. CALL LICHT(INORM,MPNORM,TYPE,IGEOMC)
  215. C
  216. C**** LICHT active les MPOVALs en *MOD
  217. C
  218. C i.e.
  219. C
  220. C SEGACT MPOVSU*MOD
  221. C SEGACT MPOVOL*MOD
  222. C SEGACT MPOVDI*MOD
  223. C SEGACT MPNORM*MOD
  224. C
  225. CALL LICHT(IDU,MPODU,TYPE,IGEOMC)
  226. C SEGACT MPODU*MOD
  227. C
  228. C USDTI initialisé a zero; USDTI = 1 / DT
  229. C
  230. NCEN=MPOVOL.VPOCHA(/1)
  231. I1=NCEN
  232. SEGINI USDTI
  233. C
  234. CALL LICHT(IRES,MPORES,TYPE,IGEOMC)
  235. CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
  236. CALL LICHT(IGN,MPGN,TYPE,IGEOMC)
  237. CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC)
  238. CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
  239. CALL LICHT(IVCO,MPOCO,TYPE,IGEOMC)
  240. CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM)
  241. CALL LICHT(ICOEV,MPCOEV,TYPE,IGEOMC)
  242. C
  243. C SEGACT MPORES*MOD
  244. C SEGACT MPRN*MOD
  245. C SEGACT MPGN*MOD
  246. C SEGACT MPRETN*MOD
  247. C SEGACT MPGAMN*MOD
  248. C SEGACT MPOCO
  249. C SEGACT MPLIM*MOD
  250. C SEGACT MPCOEV*MOD
  251. C SEGACT MPTLI*MOD
  252. C
  253. SEGINI, MPOCO1=MPOCO
  254. C
  255. C**** MPOCO1 is the "corrected cut-off" speed.
  256. C It is bigger than the given cut-off and the local
  257. C speed.
  258. C
  259. IF(IGEOMF .NE. MELEMF)THEN
  260. WRITE(IOIMP,*) 'Anomalie dans kfm13.eso'
  261. WRITE(IOIMP,*) 'Probleme de SPG'
  262. C 21 2
  263. C Données incompatibles
  264. CALL ERREUR(21)
  265. GOTO 9999
  266. ENDIF
  267. IF(IGEOMC .NE. MELEMC)THEN
  268. WRITE(IOIMP,*) 'Anomalie dans kfm13.eso'
  269. WRITE(IOIMP,*) 'Probleme de SPG'
  270. C 21 2
  271. C Données incompatibles
  272. CALL ERREUR(21)
  273. GOTO 9999
  274. ENDIF
  275. C
  276. C**** On initialise GAMMA=I+Q1.(Q2)^t
  277. C
  278. I1=5
  279. I2=NCEN
  280. SEGINI Q1
  281. SEGINI Q2
  282. C
  283. C**** On initialise le segment de travail D1
  284. C
  285. I1=5
  286. SEGINI D1
  287. SEGINI NEWLI
  288. C
  289. C**** On initialise AN0, DN0, PHI2V
  290. C
  291. I1=NCEN
  292. SEGINI AN0
  293. SEGINI DN0
  294. SEGINI PHI2V
  295. C
  296. C**** On initialise BN0
  297. C
  298. I1=5
  299. I2=NCEN
  300. SEGINI BN0
  301. C
  302. C**** On initialise BN1
  303. C
  304. I1=5
  305. I2=NCEN
  306. SEGINI BN1
  307. C
  308. C**** On initialise CN0
  309. C
  310. I1=5
  311. I2=NCEN
  312. SEGINI CN0
  313. C
  314. C**** On initialise CN1
  315. C
  316. I1=5
  317. I2=NCEN
  318. SEGINI CN1
  319. C
  320. C**** On initialise SIGMA0, SIGMAV
  321. C
  322. I1=NFAC
  323. SEGINI SIGMA0
  324. SEGINI SIGMAV
  325. C
  326. C**** RES
  327. C
  328. I1=5
  329. I2=NCEN
  330. SEGINI RES
  331. DO ICEN=1,NCEN,1
  332. DO ICOMP=1,5,1
  333. RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP)
  334. ENDDO
  335. ENDDO
  336. SEGDES MPORES
  337. C
  338. C***** DU (increment des variables conservatives)
  339. C
  340. I1=5
  341. I2=NCEN
  342. SEGINI DU
  343. C
  344. C**** Les variables conservatives
  345. C
  346. C CONS1.VAL contient les variables conservatives + gamma
  347. C i.e. rho(i) = CONS.VAL(1,i)
  348. C rhoux(i) = CONS.VAL(2,i)
  349. C rhouy(i) = CONS.VAL(3,i)
  350. C rhouz(i) = CONS.VAL(4,i)
  351. C rhoet(i) = CONS.VAL(5,i)
  352. C gam(i) = CONS.VAL(6,i)
  353. C
  354. I1=6
  355. I2=NCEN
  356. SEGINI CONS1
  357. DO ICEN=1,NCEN,1
  358. CONS1.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1)
  359. CONS1.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1)
  360. CONS1.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2)
  361. CONS1.VAL(4,ICEN)=MPGN.VPOCHA(ICEN,3)
  362. CONS1.VAL(5,ICEN)=MPRETN.VPOCHA(ICEN,1)
  363. CONS1.VAL(6,ICEN)=MPGAMN.VPOCHA(ICEN,1)
  364. ENDDO
  365. SEGDES MPRN
  366. SEGDES MPRETN
  367. SEGDES MPGN
  368. SEGDES MPGAMN
  369. C
  370. SEGINI, CONS2=CONS1
  371. C
  372. C**** BOUCLE JACOBI
  373. C
  374. DO IBLJAC=1,NJAC,1
  375. C
  376. C**** BOUCLE SUR FACEL pour le calcul de AN0, BN0
  377. C
  378. DO NLCF = 1, NFAC
  379. C
  380. C******* NLCF = numero local du centre de facel
  381. C NGCF = numero global du centre de facel
  382. C NLCF1 = numero local du centre de face
  383. C NGCEG = numero global du centre ELT "gauche"
  384. C NLCEG = numero local du centre ELT "gauche"
  385. C NGCED = numero global du centre ELT "droite"
  386. C NLCED = numero local du centre ELT "droite"
  387. C
  388. NGCEG = IPT2.NUM(1,NLCF)
  389. NGCED = IPT2.NUM(3,NLCF)
  390. NGCF = IPT2.NUM(2,NLCF)
  391. NLCF1 = MLENT2.LECT(NGCF)
  392. NLCEG = MLENT1.LECT(NGCEG)
  393. NLCED = MLENT1.LECT(NGCED)
  394. C
  395. C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
  396. C
  397. IF(NLCF .NE. NLCF1)THEN
  398. WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  399. WRITE(IOIMP,*) 'Probleme de SPG'
  400. C 21 2
  401. C Données incompatibles
  402. CALL ERREUR(21)
  403. GOTO 9999
  404. ENDIF
  405. C
  406. CNX = MPNORM.VPOCHA(NLCF,7)
  407. CNY = MPNORM.VPOCHA(NLCF,8)
  408. CNZ = MPNORM.VPOCHA(NLCF,9)
  409. CTX = MPNORM.VPOCHA(NLCF,1)
  410. CTY = MPNORM.VPOCHA(NLCF,2)
  411. CTZ = MPNORM.VPOCHA(NLCF,3)
  412. CVX = MPNORM.VPOCHA(NLCF,4)
  413. CVY = MPNORM.VPOCHA(NLCF,5)
  414. CVZ = MPNORM.VPOCHA(NLCF,6)
  415.  
  416. SURF= MPOVSU.VPOCHA(NLCF,1)
  417. C
  418. C******* Recuperation des Etats "gauche" et "droite"
  419. C
  420. C
  421. ROG = CONS2.VAL(1,NLCEG)
  422. RUXG = CONS2.VAL(2,NLCEG)
  423. RUYG = CONS2.VAL(3,NLCEG)
  424. RUZG = CONS2.VAL(4,NLCEG)
  425. RETG = CONS2.VAL(5,NLCEG)
  426. GAMG = CONS2.VAL(6,NLCEG)
  427. C
  428. UXG = RUXG / ROG
  429. UYG = RUYG / ROG
  430. UZG = RUZG / ROG
  431. C
  432. UNG = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ)
  433. UTG = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ)
  434. UVG = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ)
  435. C
  436. REC = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
  437. REC = REC / ROG
  438. PG = (GAMG - 1.0D0) * (RETG - REC)
  439. C
  440. VOLG = MPOVOL.VPOCHA(NLCEG,1)
  441. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  442. C
  443. ROD = CONS2.VAL(1,NLCED)
  444. RUXD = CONS2.VAL(2,NLCED)
  445. RUYD = CONS2.VAL(3,NLCED)
  446. RUZD = CONS2.VAL(4,NLCED)
  447. RETD = CONS2.VAL(5,NLCED)
  448. GAMD = CONS2.VAL(6,NLCED)
  449. C
  450. UXD = RUXD / ROD
  451. UYD = RUYD / ROD
  452. UZD = RUZD / ROD
  453. C
  454. UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
  455. UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
  456. UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
  457. C
  458. REC = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2)
  459. REC = REC / ROD
  460. PD = (GAMD - 1.0D0) * (RETD - REC)
  461. C
  462. VOLD = MPOVOL.VPOCHA(NLCED,1)
  463. DIAMD = MPOVDI.VPOCHA(NLCED,1)
  464.  
  465. C
  466. IF(NLCEG .EQ. NLCED)THEN
  467. C Murs ou condition limite
  468. NLLIM=MLELIM.LECT(NGCF)
  469. IF(NLLIM .EQ.0)THEN
  470. C Mur
  471. UND = -1.0D0 * UNG
  472. UTD = UTG
  473. UVD = UVG
  474.  
  475. UXD = UND * CNX + UTD * CTX + UVD * CVX
  476. UYD = UND * CNY + UTD * CTY + UVD * CVY
  477. UZD = UND * CNZ + UTD * CTZ + UVD * CVZ
  478. C
  479. RUXD = ROD * UXD
  480. RUYD = ROD * UYD
  481. RUZD = ROD * UZD
  482. C
  483. ELSE
  484. ROD = MPLIM.VPOCHA(NLLIM,1)
  485. UXD = MPLIM.VPOCHA(NLLIM,2)
  486. UYD = MPLIM.VPOCHA(NLLIM,3)
  487. UZD = MPLIM.VPOCHA(NLLIM,4)
  488. PD = MPLIM.VPOCHA(NLLIM,5)
  489. C
  490. RUXD = ROD * UXD
  491. RUYD = ROD * UYD
  492. RUZD = ROD * UZD
  493. C
  494. GAMD = GAMG
  495. C
  496. UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
  497. UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
  498. UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
  499. RETD = ((1.0D0/(GAMD - 1.0D0)) * PD)
  500. & + (0.5D0 * ROD * (UXD**2 + UYD**2 + UZD**2))
  501. VOLD = VOLG
  502. ENDIF
  503. ENDIF
  504. C
  505. C********** We check that density and pressure are positive
  506. C
  507. IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  508. & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  509. WRITE(IOIMP,*) 'FMM'
  510. C 22 0
  511. C*********** Opération malvenue. Résultat douteux
  512. C
  513. CALL ERREUR(22)
  514. IPROBL = 1
  515. GOTO 9998
  516. ENDIF
  517. C
  518. C********** Entahlpy
  519. C
  520. RHTG = RETG + PG
  521. RHTD = RETD + PD
  522. C
  523. C********** FIRST JACOBI ITERATION
  524. C
  525. IF(IBLJAC .EQ. 1)THEN
  526. C
  527. C************* We compute the interfacial state
  528. C
  529. ROF = 0.5D0 * (ROG + ROD)
  530. PF = 0.5D0 * (PG + PD)
  531. UNF = 0.5D0 * (UNG + UND)
  532. GAMF = 0.5D0 * (GAMG + GAMD)
  533. C
  534. C************* We compute SIGMA and the corrected cut-off speed
  535. C Here SIGMA is the spectral of the preconditioned Roe
  536. C matrix (|\Gamma^{-1} A| in the referred report)
  537. C
  538. UCOF = MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1))
  539. UCOF = MAX(UCOF,((UNG**2+UTG**2+UVG**2)**0.5D0))
  540. UCOF = MAX(UCOF,((UND**2+UTD**2+UVD**2)**0.5D0))
  541. IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1)) MPOCO1.VPOCHA(NLCEG
  542. $ ,1)= UCOF
  543. IF(UCOF .GT. MPOCO1.VPOCHA(NLCED,1)) MPOCO1.VPOCHA(NLCED
  544. $ ,1)= UCOF
  545. C
  546. QN = ABS(UNF)
  547. CF = (GAMF * PF / ROF)**0.5D0
  548. IF(UCOF .GE. CF)THEN
  549. PHI2 = 1.0D0
  550. ELSEIF(QN .GT. CF)THEN
  551. PHI2 = 1.0D0
  552. ELSEIF(QN .LT. UCOF)THEN
  553. PHI2 = UCOF / CF
  554. ELSE
  555. PHI2 = QN / CF
  556. ENDIF
  557. PHI2 = PHI2 * PHI2
  558. C
  559. SIGMAF = (1.0D0 - PHI2) * QN
  560. SIGMAF = SIGMAF * SIGMAF
  561. SIGMAF = SIGMAF + (4.0D0 * PHI2 * CF * CF)
  562. SIGMAF = SIGMAF**0.5D0
  563. SIGMAF = SIGMAF + ((1.0D0 + PHI2) * QN)
  564. SIGMAF = 0.5D0 * SIGMAF
  565. C
  566. C************* We store SIGMA into SIGMA0
  567. C
  568. SIGMA0.VALVEC(NLCF) = SIGMAF
  569. C
  570. C************* We compute the dual time step
  571. C
  572. UNSDT = USDTI.VALVEC(NLCEG)
  573. UNSDTF = SIGMAF / DIAMG
  574. IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG) = UNSDTF
  575. UNSDT = USDTI.VALVEC(NLCED)
  576. UNSDTF = SIGMAF / DIAMD
  577. IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCED) = UNSDTF
  578. C
  579. C************* We compute AN0
  580. C ANO in the referred report is equal to
  581. C \left(
  582. C \frac{1}{\Delta \tau_i ^0}
  583. C + \frac{1}{2{\Omega_i}}
  584. C \sum_j \sigma_{i,j}^{0} S_j
  585. C \right)
  586. C
  587. AN0.VALVEC(NLCEG) = AN0.VALVEC(NLCEG)
  588. & + (0.5D0 * SIGMAF * SURF / VOLG)
  589. IF(NLCEG .NE. NLCED)
  590. & AN0.VALVEC(NLCED) = AN0.VALVEC(NLCED)
  591. & + (0.5D0 * SIGMAF * SURF / VOLD)
  592. C
  593. C************* We compute SIGMAV
  594. C
  595. C COEF=FG
  596. COEF = (MCOORD.XCOOR(((NGCEG-1)*4)+1)
  597. & -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX
  598. COEF = COEF
  599. & + (MCOORD.XCOOR(((NGCEG-1)*4)+2)
  600. & -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY
  601. COEF = COEF
  602. & + (MCOORD.XCOOR(((NGCEG-1)*4)+3)
  603. & -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ
  604. C
  605. C COEF1=FD
  606. COEF1 = (MCOORD.XCOOR(((NGCED-1)*4)+1)
  607. & -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX
  608. COEF1 = COEF1
  609. & + (MCOORD.XCOOR(((NGCED-1)*4)+2)
  610. & -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY
  611. COEF1 = COEF1
  612. & + (MCOORD.XCOOR(((NGCED-1)*4)+3)
  613. & -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ
  614. C
  615. C COEF=|FG|+|FD|
  616. COEF = ABS(COEF)+ABS(COEF1)
  617. SIGMAF = 0.5D0 * (MPCOEV.VPOCHA(NLCEG,1)
  618. & +MPCOEV.VPOCHA(NLCED,1)) / COEF
  619. C
  620. C
  621. SIGMAV.VALVEC(NLCF) = SIGMAF
  622. C
  623. C
  624. C************* We compute the contribution of SIGMAV to
  625. C the diagonal part of the matrix to invert
  626. C
  627. C In the given report, DNO is a_i^0, i.e.
  628. C DN0 = AN0 (previously defined) +
  629. C + (\frac{1}{{\Omega_i}}
  630. C \sum_j \sigma_{V,i,j}^{0} S_j)
  631. C + 1 / \Delta t^n
  632. C
  633. COEF = SURF / VOLG
  634. DN0.VALVEC(NLCEG) = DN0.VALVEC(NLCEG)
  635. & + COEF * SIGMAF
  636. IF(NLCEG .NE. NLCED)THEN
  637. COEF = SURF / VOLD
  638. DN0.VALVEC(NLCED) = DN0.VALVEC(NLCED)
  639. & + COEF * SIGMAF
  640. ENDIF
  641. ENDIF
  642. C
  643. C********** ALL JACOBI ITERATIONS
  644. C
  645. C
  646. C********** On calcule BN1
  647. C
  648. C BN1 in the referred report is given by
  649. C \sum_j
  650. C \left\{
  651. C \frac{S_j}{2 \Omega_i} ~
  652. C \left(
  653. C \mathcal{F}(\mathbf{U}_{j,R}^{\ell})
  654. C -
  655. C 2 \sigma_{V,i,j}^{0}
  656. C \mathbf{U}_{j,R}^{\ell}
  657. C _right)
  658. C \right\}
  659. C
  660. C Central contribution
  661. C
  662. FR = (ROD * UND)
  663. FRUX = (ROD * UND * UXD) + (PD * CNX)
  664. FRUY = (ROD * UND * UYD) + (PD * CNY)
  665. FRUZ = (ROD * UND * UZD) + (PD * CNZ)
  666. FRET = (UND * RHTD)
  667. C
  668. C
  669. COEF = 0.5D0 * SURF / VOLG
  670. C
  671. BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  672. & + (COEF * FR)
  673. BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  674. & + (COEF * FRUX)
  675. BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  676. & + (COEF * FRUY)
  677. BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  678. & + (COEF * FRUZ)
  679. BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG)
  680. & + (COEF * FRET)
  681. C
  682. C Viscous contribution
  683. C
  684. COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLG
  685. C
  686. BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  687. & - (COEF * ROD)
  688. BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  689. & - (COEF * RUXD)
  690. BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  691. & - (COEF * RUYD)
  692. BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  693. & - (COEF * RUZD)
  694. BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG)
  695. & - (COEF * RETD)
  696. C
  697. C********** On calcule CN1
  698. C
  699. C CN1 in the referred report is given by
  700. C \sum_j
  701. C \left\{
  702. C \frac{S_j}{2 \Omega_i} ~
  703. C \left(
  704. C \sigma_{i,j}^{0}
  705. C \mathbf{U}_{j,R}^{\ell}
  706. C \right)
  707. C \right\}
  708. C
  709. C
  710. COEF = 0.5D0 * SURF * SIGMA0.VALVEC(NLCF) / VOLG
  711. CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG)
  712. & + (COEF * ROD)
  713. CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG)
  714. & + (COEF * RUXD)
  715. CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG)
  716. & + (COEF * RUYD)
  717. CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG)
  718. & + (COEF * RUZD)
  719. CN1.VAL(5,NLCEG) = CN1.VAL(5,NLCEG)
  720. & + (COEF * RETD)
  721. C
  722. IF(NLCEG .NE. NLCED)THEN
  723. C
  724. C********** Domain interieur
  725. C
  726. C Central contribution
  727. C
  728. FR = (ROG * UNG)
  729. FRUX = (ROG * UNG * UXG) + (PG * CNX)
  730. FRUY = (ROG * UNG * UYG) + (PG * CNY)
  731. FRUZ = (ROG * UNG * UZG) + (PG * CNZ)
  732. FRET = (UNG * RHTG)
  733. C
  734. COEF = 0.5D0 * SURF / VOLD
  735. C
  736. BN1.VAL(1,NLCED) = BN1.VAL(1,NLCED)
  737. & - (COEF * FR)
  738. BN1.VAL(2,NLCED) = BN1.VAL(2,NLCED)
  739. & - (COEF * FRUX)
  740. BN1.VAL(3,NLCED) = BN1.VAL(3,NLCED)
  741. & - (COEF * FRUY)
  742. BN1.VAL(4,NLCED) = BN1.VAL(4,NLCED)
  743. & - (COEF * FRUZ)
  744. BN1.VAL(5,NLCED) = BN1.VAL(5,NLCED)
  745. & - (COEF * FRET)
  746. C
  747. C Viscous contribution
  748. C
  749. COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLD
  750. C
  751. BN1.VAL(1,NLCED) = BN1.VAL(1,NLCED)
  752. & - (COEF * ROG)
  753. BN1.VAL(2,NLCED) = BN1.VAL(2,NLCED)
  754. & - (COEF * RUXG)
  755. BN1.VAL(3,NLCED) = BN1.VAL(3,NLCED)
  756. & - (COEF * RUYG)
  757. BN1.VAL(4,NLCED) = BN1.VAL(4,NLCED)
  758. & - (COEF * RUZG)
  759. BN1.VAL(5,NLCED) = BN1.VAL(5,NLCED)
  760. & - (COEF * RETG)
  761. C
  762. C Numerical diffusity contribution
  763. C
  764. COEF = 0.5D0 * SURF * SIGMA0.VALVEC(NLCF) / VOLD
  765. C
  766. CN1.VAL(1,NLCED) = CN1.VAL(1,NLCED)
  767. & + (COEF * ROG)
  768. CN1.VAL(2,NLCED) = CN1.VAL(2,NLCED)
  769. & + (COEF * RUXG)
  770. CN1.VAL(3,NLCED) = CN1.VAL(3,NLCED)
  771. & + (COEF * RUYG)
  772. CN1.VAL(4,NLCED) = CN1.VAL(4,NLCED)
  773. & + (COEF * RUZG)
  774. CN1.VAL(5,NLCED) = CN1.VAL(5,NLCED)
  775. & + (COEF * RETG)
  776.  
  777. ENDIF
  778. C
  779. C********** Fin boucle sur les faces
  780. C
  781. ENDDO
  782.  
  783. C
  784. C********** Loop on the centre
  785. C Computation of GAMMA
  786. C DMAT
  787. C
  788. DO ICEN=1,NCEN,1
  789. C
  790. IF(IBLJAC .EQ. 1)THEN
  791. C
  792. ROG = CONS2.VAL(1,ICEN)
  793. RUXG = CONS2.VAL(2,ICEN)
  794. RUYG = CONS2.VAL(3,ICEN)
  795. RUZG = CONS2.VAL(4,ICEN)
  796. RETG = CONS2.VAL(5,ICEN)
  797. GAMG = CONS2.VAL(6,ICEN)
  798. C
  799. UXG = RUXG / ROG
  800. UYG = RUYG / ROG
  801. UZG = RUZG / ROG
  802. C
  803. REC = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
  804. REC = REC / ROG
  805. PG = (GAMG - 1.0D0) * (RETG - REC)
  806. C
  807. C************* Note that the Cut-off is the corrected one
  808. C
  809. CG = (GAMG * PG / ROG)**0.5D0
  810. UCOG = MPOCO1.VPOCHA(ICEN,1)
  811. QG = 2.0D0 * REC / ROG
  812. QG = QG**0.5D0
  813. C
  814. IF(UCOG .GE. CG)THEN
  815. PHI2 = 1.0D0
  816. ELSEIF(QG .GT. CG)THEN
  817. PHI2 = 1.0D0
  818. ELSEIF(QG .LT. UCOG)THEN
  819. PHI2 = UCOG / CG
  820. ELSE
  821. PHI2 = QG / CG
  822. ENDIF
  823. PHI2 = PHI2 * PHI2
  824. PHI2V.VALVEC(ICEN) = PHI2
  825. C
  826. C************* We compute GAMMA
  827. C GAMMA=I+((1-PHI2)/PHI2*Q1.(Q2)^t)
  828. C
  829. C In the referred report, Q2 is the line vector
  830. C
  831. C ( \frac{q^2}{2} , -u , -v , 1)
  832. C
  833. C and Q1 is the column vector
  834. C
  835. C \frac{\gamma - 1}{c^2}
  836. C \left(
  837. C \begin{array}{c}
  838. C 1 \\
  839. C u \\
  840. C v \\
  841. C H \\
  842. C \end{array}
  843. C
  844. C
  845. Q2.VAL(1,ICEN) = (0.5D0 * QG * QG)
  846. Q2.VAL(2,ICEN) = -1.0D0 * RUXG / ROG
  847. Q2.VAL(3,ICEN) = -1.0D0 * RUYG / ROG
  848. Q2.VAL(4,ICEN) = -1.0D0 * RUZG / ROG
  849. Q2.VAL(5,ICEN) = 1.0D0
  850. C
  851. COEF1 = (GAMG - 1.0D0) / (CG*CG)
  852. COEF = (1.0D0 / COEF1) + (0.5D0 * QG * QG)
  853. C COEF=HT
  854. Q1.VAL(1,ICEN) = COEF1
  855. Q1.VAL(2,ICEN) = COEF1 * RUXG / ROG
  856. Q1.VAL(3,ICEN) = COEF1 * RUYG / ROG
  857. Q1.VAL(4,ICEN) = COEF1 * RUZG / ROG
  858. Q1.VAL(5,ICEN) = COEF1 * COEF
  859. C
  860. C************* Ebd of the computation of AN0 and DN0
  861. C
  862. AN0.VALVEC(ICEN) = AN0.VALVEC(ICEN)
  863. & + USDTI.VALVEC(ICEN) * 2.0D0 / DCFL
  864. C
  865. DN0.VALVEC(ICEN) = DN0.VALVEC(ICEN)
  866. & + (1.0D0 / DTPS)
  867. C
  868. C
  869. C************* We save BN1, CN1 into BN0 and CN0
  870. C
  871. DO ICOMP=1,5,1
  872. BN0.VAL(ICOMP,ICEN) = BN1.VAL(ICOMP,ICEN)
  873. ENDDO
  874. C
  875. DO ICOMP=1,5,1
  876. CN0.VAL(ICOMP,ICEN) = CN1.VAL(ICOMP,ICEN)
  877. ENDDO
  878. C
  879. ENDIF
  880. C
  881. C********** IBLJAC => 1
  882. C
  883. DO ICOMP=1,5,1
  884. D1.VALVEC(ICOMP) = 0.0D0
  885. ENDDO
  886. C
  887. COEF=0.0D0
  888. DO I1=1,5,1
  889. COEF = COEF
  890. & + (Q2.VAL(I1,ICEN) * (CN1.VAL(I1,ICEN)
  891. & -CN0.VAL(I1,ICEN)))
  892. ENDDO
  893. COEF = COEF
  894. & * (1.0D0 - PHI2V.VALVEC(ICEN)) / PHI2V.VALVEC(ICEN)
  895. C
  896. DO I1=1,5,1
  897. D1.VALVEC(I1) = (CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN))
  898. & + (COEF * Q1.VAL(I1,ICEN))
  899. ENDDO
  900. C
  901. C********** At this moment D1 = \Gamma \cdot (\Delta CN1)
  902. C
  903. DO ICOMP=1,5,1
  904. C
  905. C************* First increment
  906. C
  907. D1.VALVEC(ICOMP) = ( D1.VALVEC(ICOMP)
  908. & + RES.VAL(ICOMP,ICEN)
  909. & + (BN0.VAL(ICOMP,ICEN)
  910. & -BN1.VAL(ICOMP,ICEN)))
  911. & / (AN0.VALVEC(ICEN) + DN0.VALVEC(ICEN))
  912. ENDDO
  913. C
  914. C********** At this moment, in the referred report,
  915. C D1 is (\delta \mathbf{U}')^{l+1}_i in section A3
  916. C
  917. C********** Second increment
  918. C
  919. COEF=0.0D0
  920. DO I1=1,5,1
  921. COEF = COEF + (Q2.VAL(I1,ICEN) * D1.VALVEC(I1))
  922. ENDDO
  923. C
  924. C********** COEF is the scalar product between Q2 and the first
  925. C increment.
  926. C Let us multiply it by
  927. C b_i^0 / (a_i^0 + b_i^0)
  928. C
  929. COEF = COEF
  930. & * AN0.VALVEC(ICEN)
  931. & * (PHI2V.VALVEC(ICEN) - 1.0D0)
  932. COEF = COEF
  933. & / ( AN0.VALVEC(ICEN)
  934. & + (DN0.VALVEC(ICEN) * PHI2V.VALVEC(ICEN)))
  935. C
  936. DO I1=1,5,1
  937. D1.VALVEC(I1) = D1.VALVEC(I1) + COEF * Q1.VAL(I1,ICEN)
  938. ENDDO
  939. C
  940. C********** D2 = \delta \mathbf{U}
  941. C
  942. DO ICOMP=1,5,1
  943. DU.VAL(ICOMP,ICEN) = D1.VALVEC(ICOMP)
  944. CONS2.VAL(ICOMP,ICEN) = CONS1.VAL(ICOMP,ICEN)
  945. & + D1.VALVEC(ICOMP)
  946. BN1.VAL(ICOMP,ICEN) = 0.0D0
  947. CN1.VAL(ICOMP,ICEN) = 0.0D0
  948. ENDDO
  949. ENDDO
  950. C
  951. C**** Fin boucle JACOBI
  952. C
  953. ENDDO
  954. C
  955. 9998 CONTINUE
  956. C
  957. C**** Exit Data
  958. C
  959. DO ICEN=1,NCEN,1
  960. DO ICOMP=1,5,1
  961. MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
  962. ENDDO
  963. ENDDO
  964. C
  965. C
  966. C
  967. SEGSUP DU
  968. SEGSUP RES
  969. SEGSUP CONS1
  970. SEGSUP CONS2
  971. SEGSUP Q1
  972. SEGSUP Q2
  973. SEGSUP D1
  974. SEGSUP SIGMA0
  975. SEGSUP SIGMAV
  976. SEGSUP AN0
  977. SEGSUP BN0
  978. SEGSUP BN1
  979. SEGSUP CN0
  980. SEGSUP CN1
  981. SEGSUP DN0
  982. SEGSUP PHI2V
  983. SEGSUP USDTI
  984. C
  985. SEGDES MPOVSU
  986. SEGDES MPOVDI
  987. SEGDES MPOVOL
  988. SEGDES MPNORM
  989. SEGDES MPOCO
  990. SEGDES MPCOEV
  991. C
  992. SEGSUP MPOCO1
  993. C
  994. SEGSUP NEWLI
  995. C
  996. SEGDES MPODU
  997. C
  998. SEGDES IPT2
  999. SEGSUP MLENT1
  1000. SEGSUP MLENT2
  1001. C
  1002. SEGSUP MLELIM
  1003. IF(MPLIM .GT. 0) SEGDES MPLIM
  1004. C
  1005. 9999 CONTINUE
  1006. RETURN
  1007. END
  1008. C
  1009.  
  1010.  
  1011.  
  1012.  
  1013.  

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