Télécharger kfm11.eso

Retour à la liste

Numérotation des lignes :

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

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