Télécharger kfm13.eso

Retour à la liste

Numérotation des lignes :

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

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