Télécharger kfm12.eso

Retour à la liste

Numérotation des lignes :

  1. C KFM12 SOURCE BECC 06/05/15 21:15:12 5442
  2. SUBROUTINE KFM12(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
  3. & ICHPSU,ICHPDI,ICHPVO,INORM,
  4. & MELEMC,MELEMF,MELEFL,MELTFA,DTPS,DCFL,
  5. & MELLIM,ICHLIM,NJAC,ICOEF,ICOEV,
  6. & IDU,IPROBL)
  7. C************************************************************************
  8. C
  9. C PROJET : CASTEM 2000
  10. C
  11. C NOM : KFM12
  12. C
  13. C DESCRIPTION : Voir aussi KFM1
  14. C
  15. C Cas deux dimensions, gaz "calorically perfect"
  16. C Méthode de relaxation de type SGS
  17. C
  18. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  19. C
  20. C AUTEUR : T. KLOCZKO DEN/DM2S/SFME/LTMF
  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 MELEFL : 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 : 2004 : creation
  89. C
  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. C & ,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. & ,MELEMF,IDU,NFAC, NFACE
  118. & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCEG, NLCED
  119. & ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP
  120. & ,I1,I2, IBLJAC, IPROBL, IVCO, MELTFA
  121. & ,NSOUPO,IFACT,NELT, ISOUP, IELT, IFAC, IFACG, JG
  122. & ,ISOUP1, ISOUP2, ICOEF,ICOEV
  123. & , ID, IELT1, IELT2, IFAC1, IFAC2
  124. REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, RECG, RECD, PG, VOLG
  125. & , DIAMG, ROD, RUXD, RUYD, UND, RETD, GAMD, PD
  126. & , CNX, CNY, DCFL, DTPS
  127. & , SURF, SIGMAF
  128. & , UNSDT, UNSDTF, UXD, UYD
  129. & , FR, FRUX, FRUY, FRET, COEF
  130. & , CTX, CTY, RHTD, UTD, UTG, UXG, UYG, QG, CG, UCOG
  131. & , ROF, PF, UNF, GAMF, UCOF, CF, PHI2, COEF1
  132. & , QN
  133. C
  134. CHARACTER*8 TYPE
  135. LOGICAL LOGINV, LOGJAC
  136. C
  137. C**** LES INCLUDES
  138. C
  139. -INC CCOPTIO
  140. -INC SMCHPOI
  141. POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
  142. & , MPODU.MPOVAL, MPOVOL.MPOVAL
  143. & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL
  144. & , MPNORM.MPOVAL, MPLIM.MPOVAL, MPORES.MPOVAL, MPOCO.MPOVAL
  145. & , MPOCO1.MPOVAL, MELEFL.MELEME, MELEMC.MELEME
  146. & , MPCOEV.MPOVAL
  147. -INC SMCOORD
  148. -INC SMELEME
  149. -INC SMLMOTS
  150. -INC SMLENTI
  151. POINTEUR MLELIM.MLENTI, MLESUP.MLENTI
  152. C
  153. C**** POINTEUR vecteurs
  154. C
  155. SEGMENT VEC
  156. REAL*8 VALVEC(I1)
  157. ENDSEGMENT
  158. POINTEUR SIGMA0.VEC, AN0.VEC, DN0.VEC, D1.VEC
  159. & , USDTI.VEC, PHI2V.VEC, SIGMAV.VEC
  160. C
  161. C**** POINTEUR matriciels
  162. C
  163. SEGMENT MAT
  164. REAL*8 VAL(I1,I2)
  165. ENDSEGMENT
  166. POINTEUR CONS0.MAT, CONS1.MAT, BN0.MAT, BN1.MAT
  167. & , CN0.MAT, CN1.MAT, DU.MAT
  168. & , RES.MAT
  169. & , Q1.MAT, Q2.MAT
  170. C
  171. C
  172. C**** IPROBL: si different de 0, un probleme a été detecté
  173. C
  174. IPROBL=0
  175. C
  176. C**** Initialisation des MLENTI des conditions aux limites
  177. C
  178. C
  179. CALL KRIPAD(MELLIM,MLELIM)
  180. C SEGINI MLELIM
  181. C
  182. C**** Initialisation des MELEMEs
  183. C
  184. C 'CENTRE', 'FACEL'
  185. C
  186. SEGACT MELEFL
  187. NFAC = MELEFL.NUM(/2)
  188. C
  189. C**** KRIPAD pour la correspondance global/local de centre
  190. C
  191. CALL KRIPAD(MELEMC,MLENT1)
  192. C
  193. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  194. C
  195. C Si i est le numero global d'un noeud de ICEN,
  196. C MLENT1.LECT(i) contient sa position, i.e.
  197. C
  198. C I = numero global du noeud centre
  199. C MLENT1.LECT(i) = numero local du noeud centre
  200. C
  201. C MLENT1 déjà activé, i.e.
  202. C
  203. C SEGINI MLENT1
  204. C
  205. C
  206. C**** KRIPAD pour la correspondance global/local de 'FACE'
  207. C
  208. CALL KRIPAD(MELEMF,MLENT2)
  209. C SEGINI MLENT2
  210. C
  211. C
  212. C**** CHPOINTs de la table DOMAINE
  213. C
  214. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  215. CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  216. CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC)
  217. CALL LICHT(INORM,MPNORM,TYPE,IGEOMC)
  218. C
  219. C**** LICHT active les MPOVALs en *MOD
  220. C
  221. C i.e.
  222. C
  223. C SEGACT MPOVSU*MOD
  224. C SEGACT MPOVOL*MOD
  225. C SEGACT MPOVDI*MOD
  226. C SEGACT MPNORM*MOD
  227. C
  228. CALL LICHT(IDU,MPODU,TYPE,IGEOMC)
  229. C SEGACT MPODU*MOD
  230. C
  231. C USDTI initialisé a zero; USDTI = 1 / DT
  232. C
  233. NCEN=MPOVOL.VPOCHA(/1)
  234. I1=NCEN
  235. SEGINI USDTI
  236. C
  237. CALL LICHT(IRES,MPORES,TYPE,IGEOMC)
  238. CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
  239. CALL LICHT(IGN,MPGN,TYPE,IGEOMC)
  240. CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC)
  241. CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
  242. CALL LICHT(IVCO,MPOCO,TYPE,IGEOMC)
  243. CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM)
  244. CALL LICHT(ICOEV,MPCOEV,TYPE,IGEOMC)
  245. C
  246. C SEGACT MPORES*MOD
  247. C SEGACT MPRN*MOD
  248. C SEGACT MPGN*MOD
  249. C SEGACT MPRETN*MOD
  250. C SEGACT MPGAMN*MOD
  251. C SEGACT MPOCO
  252. C SEGACT MPLIM*MOD
  253. C SEGACT MPCOEV*MOD
  254. C
  255. SEGINI, MPOCO1=MPOCO
  256. C
  257. C**** MPOCO1 is the "corrected cut-off" speed.
  258. C It is bigger than the given cut-off and the local
  259. C speed.
  260. C
  261. IF(IGEOMF .NE. MELEMF)THEN
  262. WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  263. WRITE(IOIMP,*) 'Probleme de SPG'
  264. C 21 2
  265. C Données incompatibles
  266. CALL ERREUR(21)
  267. GOTO 9999
  268. ENDIF
  269. IF(IGEOMC .NE. MELEMC)THEN
  270. WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  271. WRITE(IOIMP,*) 'Probleme de SPG'
  272. C 21 2
  273. C Données incompatibles
  274. CALL ERREUR(21)
  275. GOTO 9999
  276. ENDIF
  277. C
  278. C**** On initialise GAMMA=I+(1-phi2)/phi2*Q1.(Q2)^t
  279. C
  280. I1=4
  281. I2=NCEN
  282. SEGINI Q1
  283. SEGINI Q2
  284. C
  285. C**** On initialise le segment de travail D1
  286. C
  287. I1=4
  288. SEGINI D1
  289. C
  290. C**** On initialise AN0, DN0, PHI2V
  291. C
  292. I1=NCEN
  293. SEGINI AN0
  294. SEGINI DN0
  295. SEGINI PHI2V
  296. C
  297. C**** On initialise BN0
  298. C
  299. I1=4
  300. I2=NCEN
  301. SEGINI BN0
  302. C
  303. C**** On initialise BN1
  304. C
  305. I1=4
  306. I2=NCEN
  307. SEGINI BN1
  308. C
  309. C**** On initialise CN0
  310. C
  311. I1=4
  312. I2=NCEN
  313. SEGINI CN0
  314. C
  315. C**** On initialise CN1
  316. C
  317. I1=4
  318. I2=NCEN
  319. SEGINI CN1
  320. C
  321. C**** On initialise SIGMA0, SIGMAV
  322. C
  323. I1=NFAC
  324. SEGINI SIGMA0
  325. SEGINI SIGMAV
  326. C
  327. C**** RES
  328. C
  329. I1=4
  330. I2=NCEN
  331. SEGINI RES
  332. DO ICEN=1,NCEN,1
  333. DO ICOMP=1,4,1
  334. RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP)
  335. ENDDO
  336. ENDDO
  337. SEGDES MPORES
  338. C
  339. C***** DU (increment des variables conservatives)
  340. C
  341. I1=4
  342. I2=NCEN
  343. SEGINI DU
  344. C
  345. C**** Les variables conservatives
  346. C
  347. C CONS1.VAL contient les variables conservatives + gamma
  348. C i.e. rho(i) = CONS.VAL(1,i)
  349. C rhoux(i) = CONS.VAL(2,i)
  350. C rhouy(i) = CONS.VAL(3,i)
  351. C rhoet(i) = CONS.VAL(4,i)
  352. C gam(i) = CONS.VAL(5,i)
  353. C
  354. I1=5
  355. I2=NCEN
  356. SEGINI CONS0
  357. DO ICEN=1,NCEN,1
  358. CONS0.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1)
  359. CONS0.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1)
  360. CONS0.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2)
  361. CONS0.VAL(4,ICEN)=MPRETN.VPOCHA(ICEN,1)
  362. CONS0.VAL(5,ICEN)=MPGAMN.VPOCHA(ICEN,1)
  363. ENDDO
  364. SEGDES MPRN
  365. SEGDES MPRETN
  366. SEGDES MPGN
  367. SEGDES MPGAMN
  368. C
  369. SEGINI, CONS1=CONS0
  370. C
  371. C**** Storage of l-1 values of F at the l Jacobi iteration
  372. C
  373. IPT1=MELTFA
  374. SEGACT IPT1
  375. NSOUPO= IPT1.LISOUS(/1)
  376. NSOUPO=MAX(NSOUPO,1)
  377. JG=NSOUPO
  378. SEGINI MLESUP
  379. IFACT=0
  380. IF(NSOUPO .EQ. 1)THEN
  381. MLESUP.LECT(1)=IPT1
  382. NFACE=IPT1.NUM(/1)
  383. NELT=IPT1.NUM(/2)
  384. IFACT=IFACT+(NFACE*NELT)
  385. ELSE
  386. DO ISOUP = 1, NSOUPO, 1
  387. IPT2=IPT1.LISOUS(ISOUP)
  388. MLESUP.LECT(ISOUP)=IPT2
  389. SEGACT IPT2
  390. NFACE=IPT2.NUM(/1)
  391. NELT=IPT2.NUM(/2)
  392. IFACT=IFACT+(NFACE*NELT)
  393. ENDDO
  394. ENDIF
  395. C
  396. C***********************************************************************
  397. C**** FIRST JACOBI ITERATION *******************************************
  398. C***********************************************************************
  399. C
  400. C
  401. C**** LOOP on ELTFA to compute:
  402. C
  403. C AN0_i = 1/(2 \Omega_i) \sum_j \sigma_i,j S_j
  404. C SIGMA0_j
  405. C The cut-off in each cell
  406. C The reciprocal of the dual time step USDTI
  407. C
  408. ICEN=0
  409. IFACG=0
  410. DO ISOUP=1,NSOUPO,1
  411. IPT2=MLESUP.LECT(ISOUP)
  412. NFACE=IPT2.NUM(/1)
  413. NELT=IPT2.NUM(/2)
  414. DO IELT=1,NELT,1
  415. ICEN=ICEN+1
  416. DO IFAC=1,NFACE,1
  417. C IFACG increase with IFAC,IELT,ISOUP
  418. IFACG = IFACG+1
  419. NGCF = IPT2.NUM(IFAC,IELT)
  420. NLCF = MLENT2.LECT(NGCF)
  421. C
  422. C************* NLCF = numero local du centre de facel
  423. C NGCF = numero global du centre de facel
  424. C NGCEG = numero global du centre ELT "gauche"
  425. C NLCEG = numero local du centre ELT "gauche"
  426. C NGCED = numero global du centre ELT "droite"
  427. C NLCED = numero local du centre ELT "droite"
  428. C
  429. NGCEG = MELEFL.NUM(1,NLCF)
  430. NGCED = MELEFL.NUM(3,NLCF)
  431. NLCEG = MLENT1.LECT(NGCEG)
  432. NLCED = MLENT1.LECT(NGCED)
  433. C
  434. CNX = MPNORM.VPOCHA(NLCF,1)
  435. CNY = MPNORM.VPOCHA(NLCF,2)
  436. CTX = -1.0D0 * CNY
  437. CTY = CNX
  438. C
  439. C************* Left or right?
  440. C
  441. IF(NLCEG .EQ. ICEN)THEN
  442. LOGINV=.FALSE.
  443. ELSEIF(NLCED .EQ. ICEN)THEN
  444. LOGINV=.TRUE.
  445. ELSE
  446. WRITE(IOIMP,*) 'Anomalie dans kfm13.eso'
  447. WRITE(IOIMP,*) 'Probleme de SPG'
  448. C 21 2
  449. C Données incompatibles
  450. CALL ERREUR(21)
  451. GOTO 9999
  452. ENDIF
  453. C************* If right, inversion
  454. IF(LOGINV) THEN
  455. NLCED=NLCEG
  456. NLCEG=ICEN
  457. CNX = -1.0D0*CNX
  458. CNY = -1.0D0*CNY
  459. CTX = -1.0D0*CTX
  460. CTY = -1.0D0*CTY
  461. ENDIF
  462. C
  463. SURF = MPOVSU.VPOCHA(NLCF,1)
  464. C
  465. C************* Récupération des Etats "gauche" et "droite"
  466. C
  467. C
  468. ROG = CONS0.VAL(1,NLCEG)
  469. RUXG = CONS0.VAL(2,NLCEG)
  470. RUYG = CONS0.VAL(3,NLCEG)
  471. UXG = RUXG / ROG
  472. UYG = RUYG / ROG
  473. UNG = (UXG * CNX) + (UYG * CNY)
  474. UTG = (UXG * CTX) + (UYG * CTY)
  475. RETG = CONS0.VAL(4,NLCEG)
  476. GAMG = CONS0.VAL(5,NLCEG)
  477. RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
  478. RECG = RECG / ROG
  479. PG = (GAMG - 1.0D0) * (RETG - RECG)
  480. VOLG = MPOVOL.VPOCHA(NLCEG,1)
  481. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  482. C
  483. ROD = CONS0.VAL(1,NLCED)
  484. RUXD = CONS0.VAL(2,NLCED)
  485. RUYD = CONS0.VAL(3,NLCED)
  486. UXD = RUXD / ROD
  487. UYD = RUYD / ROD
  488. RETD = CONS0.VAL(4,NLCED)
  489. GAMD = CONS0.VAL(5,NLCED)
  490. RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
  491. RECD = RECD / ROD
  492. PD = (GAMD - 1.0D0) * (RETD - RECD)
  493. UND = (UXD * CNX) + (UYD * CNY)
  494. UTD = (UXD * CTX) + (UYD * CTY)
  495. C
  496. IF(NLCEG .EQ. NLCED)THEN
  497. C Murs ou condition limite
  498. NLLIM=MLELIM.LECT(NGCF)
  499. IF(NLLIM .EQ.0)THEN
  500. C Mur
  501. UND = -1.0D0 * UNG
  502. UTD = UTG
  503. UXD = UND * CNX + UTD * CTX
  504. UYD = UND * CNY + UTD * CTY
  505. RUXD = ROD*UXD
  506. RUYD = ROD*UYD
  507. ELSE
  508. ROD = MPLIM.VPOCHA(NLLIM,1)
  509. UXD = MPLIM.VPOCHA(NLLIM,2)
  510. UYD = MPLIM.VPOCHA(NLLIM,3)
  511. RUXD = ROD*UXD
  512. RUYD = ROD*UYD
  513. PD = MPLIM.VPOCHA(NLLIM,4)
  514. GAMD = GAMG
  515. UND = (UXD * CNX) + (UYD * CNY)
  516. UTD = (UXD * CTX) + (UYD * CTY)
  517. RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
  518. & +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
  519. ENDIF
  520. ENDIF
  521. C
  522. C************* We check that density and pressure are positive
  523. C
  524. IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  525. & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  526. WRITE(IOIMP,*) 'FMM'
  527. WRITE(IOIMP,*) 'TEST'
  528. C 22 0
  529. C**************** Opération malvenue. Résultat douteux
  530. C
  531. CALL ERREUR(22)
  532. IPROBL = 1
  533. GOTO 9998
  534. ENDIF
  535. C
  536. C************* We compute the interfacial state
  537. C
  538. ROF = 0.5D0*(ROG+ROD)
  539. PF = 0.5D0*(PG+PD)
  540. UNF = 0.5D0*(UNG+UND)
  541. GAMF = 0.5D0*(GAMG+GAMD)
  542. C
  543. C************* We compute SIGMA and the corrected cut-off speed
  544. C
  545. UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1))
  546. UCOF=MAX(UCOF,(((UNG*UNG)+(UTG*UTG))**0.5D0))
  547. UCOF=MAX(UCOF,(((UND*UND)+(UTD*UTD))**0.5D0))
  548. IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1))
  549. & MPOCO1.VPOCHA(NLCEG,1)= UCOF
  550. C
  551. QN=ABS(UNF)
  552. CF=(GAMF*PF/ROF)**0.5D0
  553. IF(UCOF .GE. CF)THEN
  554. PHI2 = 1.0D0
  555. ELSEIF(QN .GT. CF)THEN
  556. PHI2 = 1.0D0
  557. ELSEIF(QN .LT. UCOF)THEN
  558. PHI2 = UCOF / CF
  559. ELSE
  560. PHI2 = QN / CF
  561. ENDIF
  562. PHI2=PHI2*PHI2
  563. C
  564. SIGMAF=(1.0D0-PHI2)*QN
  565. SIGMAF=SIGMAF*SIGMAF
  566. SIGMAF=SIGMAF+(4.0D0*PHI2*CF*CF)
  567. SIGMAF=SIGMAF**0.5D0
  568. SIGMAF=SIGMAF+((1.0D0+PHI2)*QN)
  569. SIGMAF=0.5D0*SIGMAF
  570. C
  571. SIGMA0.VALVEC(NLCF)=SIGMAF
  572. C
  573. C************* We compute AN0
  574. C
  575. AN0.VALVEC(NLCEG)=AN0.VALVEC(NLCEG)
  576. & +(0.5D0*SIGMAF*SURF/VOLG)
  577. C
  578. C************* We compute the dual time step
  579. C
  580. UNSDT=USDTI.VALVEC(NLCEG)
  581. UNSDTF = SIGMAF / DIAMG
  582. IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG)=UNSDTF
  583. C
  584. C************* We compute SIGMAV
  585. C
  586. C COEF=FG
  587. COEF=(MCOORD.XCOOR(((NGCEG-1)*3)+1)-
  588. & MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX
  589. COEF=COEF+
  590. & ((MCOORD.XCOOR(((NGCEG-1)*3)+2)-
  591. & MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY)
  592. C COEF1=FD
  593. COEF1=(MCOORD.XCOOR(((NGCED-1)*3)+1)-
  594. & MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX
  595. COEF1=COEF1+
  596. & ((MCOORD.XCOOR(((NGCED-1)*3)+2)-
  597. & MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY)
  598. C COEF=|FG|+|FD|
  599. COEF=ABS(COEF)+ABS(COEF1)
  600. SIGMAF=0.5D0*(MPCOEV.VPOCHA(NLCEG,1)+
  601. & MPCOEV.VPOCHA(NLCED,1))/COEF
  602. SIGMAV.VALVEC(NLCF)=SIGMAF
  603. C
  604. C************* We compute the contribution of SIGMAV to
  605. C the diagonal part of the matrix to invert
  606. C
  607. C In the given report, DNO is a_i^0, i.e.
  608. C DN0 = AN0 (previously defined) +
  609. C + (\frac{1}{{\Omega_i}}
  610. C \sum_j \sigma_{V,i,j}^{0} S_j)
  611. C + 1 / \Delta t^n
  612. C
  613. COEF=SURF/VOLG
  614. DN0.VALVEC(NLCEG)=DN0.VALVEC(NLCEG)+COEF*SIGMAF
  615. C
  616. C
  617. C************* We compute the flux
  618. C FLU.VALVEC(1:4,IFAC) flux in IFAC face
  619. C 1 -> mass
  620. C 2 -> momentum (x)
  621. C 3 -> momentum (y)
  622. C 4 -> energy
  623. C
  624. RHTD = RETD+PD
  625. FR = (ROD*UND)
  626. FRUX = (ROD*UND*UXD)+(PD*CNX)
  627. FRUY = (ROD*UND*UYD)+(PD*CNY)
  628. FRET = (UND*RHTD)
  629. C
  630. COEF = SURF/(2.0D0*VOLG)
  631. C
  632. BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
  633. & + FR * COEF
  634. BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
  635. & + FRUX * COEF
  636. BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
  637. & + FRUY * COEF
  638. BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
  639. & + FRET * COEF
  640. C
  641. C************* Viscous contribution
  642. C
  643. COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
  644. C
  645. BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
  646. & - (COEF * ROD)
  647. BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
  648. & - (COEF * RUXD)
  649. BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
  650. & - (COEF * RUYD)
  651. BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
  652. & - (COEF * RETD)
  653. C
  654. C************* On calcule CN0
  655. C
  656. C CN0 in the referred report is given by
  657. C \sum_j
  658. C \left\{
  659. C \frac{S_j}{2 \Omega_i} ~
  660. C \left(
  661. C \sigma_{i,j}^{0}
  662. C \mathbf{U}_{j,R}^{n+1,m}
  663. C \right)
  664. C \right\}
  665. C
  666. C
  667. COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
  668. C
  669. CN0.VAL(1,NLCEG) = CN0.VAL(1,NLCEG)
  670. & + (COEF * ROD)
  671. CN0.VAL(2,NLCEG) = CN0.VAL(2,NLCEG)
  672. & + (COEF * RUXD)
  673. CN0.VAL(3,NLCEG) = CN0.VAL(3,NLCEG)
  674. & + (COEF * RUYD)
  675. CN0.VAL(4,NLCEG) = CN0.VAL(4,NLCEG)
  676. & + (COEF * RETD)
  677.  
  678. ENDDO
  679. C
  680. C********** End of loop on the faces of the element
  681. C
  682. C N.B.: ICEN = NLCEG
  683. C
  684. C
  685. C Computation of the low-Mach preconditioning Matrix
  686. C
  687. CG=(GAMG*PG/ROG)**0.5D0
  688. UCOG=MPOCO1.VPOCHA(ICEN,1)
  689. QG=2.0D0*RECG/ROG
  690. QG=QG**0.5D0
  691. C
  692. IF(UCOG .GE. CG)THEN
  693. PHI2 = 1.0D0
  694. ELSEIF(QG .GT. CG)THEN
  695. PHI2 = 1.0D0
  696. ELSEIF(QG .LT. UCOG)THEN
  697. PHI2 = UCOG / CG
  698. ELSE
  699. PHI2 = QG / CG
  700. ENDIF
  701. PHI2=PHI2*PHI2
  702. PHI2V.VALVEC(ICEN)=PHI2
  703. C
  704. Q2.VAL(1,ICEN)=(0.5D0*QG*QG)
  705. Q2.VAL(2,ICEN)=-1.0D0*RUXG/ROG
  706. Q2.VAL(3,ICEN)=-1.0D0*RUYG/ROG
  707. Q2.VAL(4,ICEN)=1.0D0
  708. C
  709. COEF1=(GAMG-1.0D0)/(CG*CG)
  710. COEF=(1.0D0/COEF1)+(0.5D0*QG*QG)
  711. C
  712. Q1.VAL(1,ICEN)=COEF1
  713. Q1.VAL(2,ICEN)=COEF1*RUXG/ROG
  714. Q1.VAL(3,ICEN)=COEF1*RUYG/ROG
  715. Q1.VAL(4,ICEN)=COEF1*COEF
  716. C
  717. C
  718. C Computation of the diagonal coefficients
  719. C
  720. AN0.VALVEC(ICEN)=AN0.VALVEC(ICEN)
  721. C & +(USDTI.VALVEC(ICEN)*2.0D0/DCFL)
  722. C
  723. DN0.VALVEC(ICEN)=DN0.VALVEC(ICEN)+(1.0D0/DTPS)
  724. & +(USDTI.VALVEC(ICEN)*2.0D0/DCFL)
  725. C
  726. ENDDO
  727. ENDDO
  728. C
  729. C***********************************************************************
  730. C**** END FIRST JACOBI ITERATION ***************************************
  731. C***********************************************************************
  732. C
  733. C***********************************************************************
  734. C**** JACOBI LOOP ******************************************************
  735. C***********************************************************************
  736. C
  737. C
  738. LOGJAC=.TRUE.
  739. C
  740. DO IBLJAC=1,NJAC,1
  741. C
  742. C
  743. C Si ICOEF=2 (i.e. si on a "LJACOF" (voir kfm1.eso))
  744. C alors relaxation SGS avec uniquement des balayages "forward"
  745. C
  746. C Si ICOEF=3 (i.e. si on a "LJACOB")
  747. C alors relaxation SGS avec uniquement des balayages "backward"
  748. C
  749. C Sinon (i.e. si on a "LJACOFB")
  750. C alors relaxation SGS avec des balayages "forward" et "backward"
  751. C
  752. IF(ICOEF .EQ. 2)THEN
  753. LOGJAC=.TRUE.
  754. ELSEIF(ICOEF .EQ. 3)THEN
  755. LOGJAC=.FALSE.
  756. ELSE
  757. LOGINV=LOGJAC
  758. LOGJAC=(IBLJAC / 2 * 2) .EQ. IBLJAC
  759. IF(LOGJAC)THEN
  760. LOGJAC= .NOT. LOGINV
  761. ELSE
  762. LOGJAC= LOGINV
  763. ENDIF
  764. ENDIF
  765. C
  766. IF(LOGJAC)THEN
  767. ICEN=0
  768. IFACG=0
  769. ID=1
  770. ISOUP1=1
  771. ISOUP2=NSOUPO
  772. ELSE
  773. ICEN=NCEN+1
  774. IFACG=IFACT+1
  775. ID=-1
  776. ISOUP1=NSOUPO
  777. ISOUP2=1
  778. ENDIF
  779. C
  780. DO ISOUP=ISOUP1,ISOUP2,ID
  781. IPT2=MLESUP.LECT(ISOUP)
  782. NFACE=IPT2.NUM(/1)
  783. NELT=IPT2.NUM(/2)
  784. C
  785. IF(LOGJAC)THEN
  786. IELT1=1
  787. IELT2=NELT
  788. IFAC1=1
  789. IFAC2=NFACE
  790. ELSE
  791. IELT2=1
  792. IELT1=NELT
  793. IFAC2=1
  794. IFAC1=NFACE
  795. ENDIF
  796. C
  797. DO IELT=IELT1,IELT2,ID
  798. ICEN=ICEN+ID
  799. DO IFAC=IFAC1,IFAC2,ID
  800. C IFACG increase with IFAC,IELT,ISOUP
  801. IFACG = IFACG+ID
  802. NGCF = IPT2.NUM(IFAC,IELT)
  803. NLCF = MLENT2.LECT(NGCF)
  804. NGCEG = MELEFL.NUM(1,NLCF)
  805. NGCED = MELEFL.NUM(3,NLCF)
  806. NLCEG = MLENT1.LECT(NGCEG)
  807. NLCED = MLENT1.LECT(NGCED)
  808. C
  809. CNX = MPNORM.VPOCHA(NLCF,1)
  810. CNY = MPNORM.VPOCHA(NLCF,2)
  811. CTX = -1.0D0 * CNY
  812. CTY = CNX
  813. C
  814. C**************** Left or right?
  815. C
  816. IF(NLCEG .EQ. ICEN)THEN
  817. LOGINV=.FALSE.
  818. ELSEIF(NLCED .EQ. ICEN)THEN
  819. LOGINV=.TRUE.
  820. ELSE
  821. WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
  822. WRITE(IOIMP,*) 'Probleme de SPG'
  823. C 21 2
  824. C Données incompatibles
  825. CALL ERREUR(21)
  826. GOTO 9999
  827. ENDIF
  828. C**************** If right, inversion
  829. IF(LOGINV) THEN
  830. NLCED=NLCEG
  831. NLCEG=ICEN
  832. CNX = -1.0D0*CNX
  833. CNY = -1.0D0*CNY
  834. CTX = -1.0D0*CTX
  835. CTY = -1.0D0*CTY
  836. ENDIF
  837. C
  838. SURF = MPOVSU.VPOCHA(NLCF,1)
  839. SIGMAF = SIGMA0.VALVEC(NLCF)
  840. C
  841. C**************** Recuperation des Etats "gauche" et "droite"
  842. C
  843. ROG = CONS1.VAL(1,NLCEG)
  844. RUXG = CONS1.VAL(2,NLCEG)
  845. RUYG = CONS1.VAL(3,NLCEG)
  846. ROG = CONS1.VAL(1,NLCEG)
  847. RUXG = CONS1.VAL(2,NLCEG)
  848. RUYG = CONS1.VAL(3,NLCEG)
  849. UXG = RUXG / ROG
  850. UYG = RUYG / ROG
  851. UNG = (UXG * CNX) + (UYG * CNY)
  852. UTG = (UXG * CTX) + (UYG * CTY)
  853. RETG = CONS1.VAL(4,NLCEG)
  854. GAMG = CONS1.VAL(5,NLCEG)
  855. RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
  856. RECG = RECG / ROG
  857. PG = (GAMG - 1.0D0) * (RETG - RECG)
  858. VOLG = MPOVOL.VPOCHA(NLCEG,1)
  859. C
  860. ROD = CONS1.VAL(1,NLCED)
  861. RUXD = CONS1.VAL(2,NLCED)
  862. RUYD = CONS1.VAL(3,NLCED)
  863. UXD = RUXD / ROD
  864. UYD = RUYD / ROD
  865. RETD = CONS1.VAL(4,NLCED)
  866. GAMD = CONS1.VAL(5,NLCED)
  867. RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
  868. RECD = RECD / ROD
  869. PD = (GAMD - 1.0D0) * (RETD - RECD)
  870. UND = (UXD * CNX) + (UYD * CNY)
  871. UTD = (UXD * CTX) + (UYD * CTY)
  872. C
  873. IF(NLCEG .EQ. NLCED)THEN
  874. C Murs ou condition limite
  875. NLLIM=MLELIM.LECT(NGCF)
  876. IF(NLLIM .EQ.0)THEN
  877. C Mur
  878. UND = -1.0D0 * UNG
  879. UTD = UTG
  880. UXD = UND * CNX + UTD * CTX
  881. UYD = UND * CNY + UTD * CTY
  882. RUXD = ROD*UXD
  883. RUYD = ROD*UYD
  884. ELSE
  885. ROD = MPLIM.VPOCHA(NLLIM,1)
  886. UXD = MPLIM.VPOCHA(NLLIM,2)
  887. UYD = MPLIM.VPOCHA(NLLIM,3)
  888. RUXD = ROD*UXD
  889. RUYD = ROD*UYD
  890. PD = MPLIM.VPOCHA(NLLIM,4)
  891. GAMD = GAMG
  892. UND = (UXD * CNX) + (UYD * CNY)
  893. UTD = (UXD * CTX) + (UYD * CTY)
  894. RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
  895. & +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
  896. ENDIF
  897. ENDIF
  898. C
  899. C**************** We check that density and pressure are positive
  900. C
  901. IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  902. & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  903. WRITE(IOIMP,*) 'FMM'
  904. WRITE(IOIMP,*) 'ITERATION JAC',IBLJAC
  905. WRITE(IOIMP,*) 'ITERATION ICEN',ICEN
  906. WRITE(IOIMP,*) ROG,PG
  907. WRITE(IOIMP,*) ROD,PD
  908. C
  909. C 22 0
  910. C******************* Opération malvenue. Résultat douteux
  911. C
  912. CALL ERREUR(22)
  913. IPROBL = 1
  914. GOTO 9998
  915. ENDIF
  916. C
  917. C
  918. C**************** We compute the flux
  919. C FLU.VALVEC(1:4,IFAC) flux in IFAC face
  920. C 1 -> mass
  921. C 2 -> momentum (x)
  922. C 3 -> momentum (y)
  923. C 4 -> energy
  924. C
  925. RHTD = RETD+PD
  926. FR = (ROD*UND)
  927. FRUX = (ROD*UND*UXD)+(PD*CNX)
  928. FRUY = (ROD*UND*UYD)+(PD*CNY)
  929. FRET = (UND*RHTD)
  930. C
  931. COEF = SURF/(2.0D0*VOLG)
  932. C
  933. BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  934. & + FR * COEF
  935. BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  936. & + FRUX * COEF
  937. BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  938. & + FRUY * COEF
  939. BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  940. & + FRET * COEF
  941. C
  942. C**************** Viscous contribution
  943. C
  944. COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
  945. C
  946. BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  947. & - (COEF * ROD)
  948. BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  949. & - (COEF * RUXD)
  950. BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  951. & - (COEF * RUYD)
  952. BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  953. & - (COEF * RETD)
  954. C
  955. C**************** On calcule DIS1
  956. C
  957. C DIS1 in the referred report is given by
  958. C \sum_j
  959. C \left\{
  960. C \frac{S_j}{2 \Omega_i} ~
  961. C \left(
  962. C \sigma_{i,j}^{0}
  963. C \mathbf{U}_{j,R}^{n+1,m,l}
  964. C \right)
  965. C \right\}
  966. C
  967. C
  968. COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
  969. C
  970. CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG)
  971. & + (COEF * ROD)
  972. CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG)
  973. & + (COEF * RUXD)
  974. CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG)
  975. & + (COEF * RUYD)
  976. CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG)
  977. & + (COEF * RETD)
  978.  
  979. ENDDO
  980. C
  981. C************* End of the loop on the faces of the element
  982. C
  983. DO ICOMP=1,4,1
  984. D1.VALVEC(ICOMP) = 0.0D0
  985. ENDDO
  986. C
  987. COEF=0.0D0
  988. DO I1=1,4,1
  989. COEF=COEF+(Q2.VAL(I1,ICEN)
  990. & *(CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN)))
  991. ENDDO
  992. COEF=COEF*(1.0D0-PHI2V.VALVEC(ICEN))
  993. & /PHI2V.VALVEC(ICEN)
  994. C
  995. DO I1=1,4,1
  996. D1.VALVEC(I1)=(CN1.VAL(I1,ICEN)
  997. & - CN0.VAL(I1,ICEN))
  998. & +(COEF*Q1.VAL(I1,ICEN))
  999. c write(6,*)D1.VALVEC(I1)
  1000. ENDDO
  1001. C
  1002. C************* At this moment D1 = \Gamma \cdot (\Delta CN1)
  1003. C
  1004. DO ICOMP=1,4,1
  1005. C
  1006. C************* First increment
  1007. C
  1008. D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP)
  1009. & +RES.VAL(ICOMP,ICEN)
  1010. & +(BN0.VAL(ICOMP,ICEN)-BN1.VAL(ICOMP,ICEN)))
  1011. & /(AN0.VALVEC(ICEN)+DN0.VALVEC(ICEN))
  1012. ENDDO
  1013. C
  1014. C************* At this moment, in the referred report,
  1015. C D1 is (\delta \mathbf{U}')^{l+1}_i in section A3
  1016. C
  1017. C************* Second increment
  1018. C
  1019. COEF=0.0D0
  1020. DO I1=1,4,1
  1021. COEF=COEF+(Q2.VAL(I1,ICEN)*D1.VALVEC(I1))
  1022. ENDDO
  1023. C
  1024. C************* COEF is the scalar product between Q2 and the first
  1025. C increment.
  1026. C Let us multiply it by
  1027. C b_i^0 / (a_i^0 + b_i^0)
  1028. C
  1029. COEF=COEF*AN0.VALVEC(ICEN)
  1030. & *(PHI2V.VALVEC(ICEN)-1.0D0)
  1031. COEF=COEF/(AN0.VALVEC(ICEN)+(DN0.VALVEC(ICEN)
  1032. & *PHI2V.VALVEC(ICEN)))
  1033. C
  1034. DO I1=1,4,1
  1035. D1.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN)
  1036. ENDDO
  1037. C
  1038. C************* D1 = \delta \mathbf{U}
  1039. C
  1040. DO ICOMP=1,4,1
  1041. DU.VAL(ICOMP,ICEN) = D1.VALVEC(ICOMP)
  1042. C
  1043. CONS1.VAL(ICOMP,ICEN) = CONS0.VAL(ICOMP,ICEN)
  1044. & + D1.VALVEC(ICOMP)
  1045. C
  1046. BN1.VAL(ICOMP,ICEN) = 0.0D0
  1047. CN1.VAL(ICOMP,ICEN) = 0.0D0
  1048. ENDDO
  1049. C
  1050. C********** End of the loop on the elements
  1051. ENDDO
  1052. C
  1053. C******* End of the loop on the Domains
  1054. ENDDO
  1055. C
  1056. C**** End of the Jacobi loop
  1057. ENDDO
  1058. C
  1059. 9998 CONTINUE
  1060. DO ICEN=1,NCEN,1
  1061. DO ICOMP=1,4,1
  1062. MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
  1063. ENDDO
  1064. ENDDO
  1065. C
  1066. SEGSUP DU
  1067. SEGSUP RES
  1068. SEGSUP CONS0
  1069. SEGSUP CONS1
  1070. SEGSUP Q1
  1071. SEGSUP Q2
  1072. SEGSUP SIGMA0
  1073. SEGSUP AN0
  1074. SEGSUP AN0
  1075. SEGSUP BN0
  1076. SEGSUP BN1
  1077. SEGSUP CN0
  1078. SEGSUP CN1
  1079. SEGSUP DN0
  1080. SEGSUP PHI2V
  1081. SEGSUP USDTI
  1082. C
  1083. SEGDES MPOVSU
  1084. SEGDES MPOVDI
  1085. SEGDES MPOVOL
  1086. SEGDES MPNORM
  1087. SEGDES MPOCO
  1088. SEGDES MPCOEV
  1089. C
  1090. SEGSUP MPOCO1
  1091. C
  1092. SEGDES MPODU
  1093. C
  1094. SEGDES MELEFL
  1095. SEGSUP MLENT1
  1096. SEGSUP MLENT2
  1097. C
  1098. SEGSUP MLELIM
  1099. IF(MPLIM .GT. 0) SEGDES MPLIM
  1100. C
  1101. DO ISOUP=1,NSOUPO,1
  1102. IPT2=MLESUP.LECT(ISOUP)
  1103. SEGDES IPT2
  1104. ENDDO
  1105. IPT2 = MELTFA
  1106. IF(ISOUP .GT. 1) SEGDES IPT2
  1107. SEGSUP MLESUP
  1108. C
  1109. 9999 CONTINUE
  1110. RETURN
  1111. END
  1112. C
  1113.  
  1114.  
  1115.  
  1116.  

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