Télécharger kfm11.eso

Retour à la liste

Numérotation des lignes :

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

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