Télécharger kfm14.eso

Retour à la liste

Numérotation des lignes :

  1. C KFM14 SOURCE BECC 06/05/15 21:15:28 5442
  2. C KFM14 SOURCE CHAT 05/01/13 00:55:31 5004
  3. SUBROUTINE KFM14(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
  4. & ICHPSU,ICHPDI,ICHPVO,INORM,
  5. & MELEMC,MELEMF,MELEFL,MELTFA,DTPS,DCFL,
  6. & MELLIM,ICHLIM,NJAC,ICOEF,ICOEV,
  7. & IDU,IPROBL)
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : KFM14
  13. C
  14. C DESCRIPTION : Voir aussi KFM1 et KFM12
  15. C
  16. C Cas 3D, gaz "calorically perfect"
  17. C relaxation de type SGS
  18. C
  19. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  20. C
  21. C AUTEUR : T. KLOCZKO, DEN/DM2S/SFME/LTMF
  22. C
  23. C************************************************************************
  24. C ENTREES
  25. C
  26. C
  27. C 1) Pointeurs des CHPOINTs
  28. C
  29. C IRES : CHPOINT 'CENTRE' contenant le residu
  30. C
  31. C IRN : CHPOINT 'CENTRE' contenant la masse volumique
  32. C
  33. C IGN : CHPOINT 'CENTRE' contenant la q.d.m.
  34. C
  35. C IRETN : CHPOINT 'CENTRE' contenant l'energie totale
  36. C
  37. C IGAMN : CHPOINT 'CENTRE' contenant le gamma
  38. C
  39. C IVCO : CHPOINT 'CENTRE' contenant le cut-off
  40. C
  41. C ICHLIM : CHPOINT conditions aux bords
  42. C
  43. C ICOEV : CHPOINT 'CENTRE' contenant le (gamma - 1) K / (R \rho)
  44. C coeff pour calculer le ray. spect. de la matrice
  45. C visc.
  46. C
  47. C 3) Pointeurs de CHPOINTs de la table DOMAINE
  48. C
  49. C ICHPSU : CHPOINT "FACE" contenant la surface des faces
  50. C
  51. C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
  52. C de chaque element
  53. C
  54. C ICHPVO : CHPOINT "CENTRE" contenant le volume
  55. C de chaque element
  56. C
  57. C INORM : CHPOINT "CENTRE" contenant le normales aux faces
  58. C
  59. C
  60. C 4) Pointeurs de MELEME de la table DOMAINE
  61. C
  62. C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
  63. C
  64. C MELEMF : MELEME 'FACE' du SPG des FACES
  65. C
  66. C MELEFL : MELEME 'FACEL' du connectivité FACES -> ELEM
  67. C
  68. C MELLIM : MELEME SPG conditions aux bords
  69. C
  70. C 5)
  71. C
  72. C DCFL = le double de la CFL
  73. C
  74. C DTPS = pas de temps physique
  75. C
  76. C NJAC = nombre d'iterations dans le jacobien
  77. C
  78. C
  79. C SORTIES
  80. C
  81. C IDU : resultat
  82. C
  83. C IPROBL : si > 0, il y a un probleme
  84. C
  85. C************************************************************************
  86. C
  87. C HISTORIQUE (Anomalies et modifications éventuelles)
  88. C
  89. C HISTORIQUE : 2004 création
  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, RUZG, RETG, GAMG, RHTG, RECG
  125. & , ROD, RUXD, RUYD, RUZD, RETD, GAMD, RHTD, RECD
  126. & , UXG, UYG, UZG, UNG, UTG, UVG, PG, VOLG, DIAMG, SURF
  127. & , UXD, UYD, UZD, UND, UTD, UVD, PD
  128. & , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ
  129. & , DCFL, DTPS, UNSDT, UNSDTF
  130. & , FR, FRUX, FRUY, FRUZ, FRET
  131. & , QG, CG, UCOG, ROF, PF, UNF, GAMF, UCOF, CF, SIGMAF
  132. & , QN, PHI2, COEF, COEF1
  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 dans kfm14.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 dans kfm14.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=5
  281. I2=NCEN
  282. SEGINI Q1
  283. SEGINI Q2
  284. C
  285. C**** On initialise le segment de travail D1
  286. C
  287. I1=5
  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=5
  300. I2=NCEN
  301. SEGINI BN0
  302. C
  303. C**** On initialise BN1
  304. C
  305. I1=5
  306. I2=NCEN
  307. SEGINI BN1
  308. C
  309. C**** On initialise CN0
  310. C
  311. I1=5
  312. I2=NCEN
  313. SEGINI CN0
  314. C
  315. C**** On initialise CN1
  316. C
  317. I1=5
  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=5
  330. I2=NCEN
  331. SEGINI RES
  332. DO ICEN=1,NCEN,1
  333. DO ICOMP=1,5,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=5
  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 rhouz(i) = CONS.VAL(4,i)
  352. C rhoet(i) = CONS.VAL(5,i)
  353. C gam(i) = CONS.VAL(6,i)
  354. C
  355. I1=6
  356. I2=NCEN
  357. SEGINI CONS0
  358. DO ICEN=1,NCEN,1
  359. CONS0.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1)
  360. CONS0.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1)
  361. CONS0.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2)
  362. CONS0.VAL(4,ICEN)=MPGN.VPOCHA(ICEN,3)
  363. CONS0.VAL(5,ICEN)=MPRETN.VPOCHA(ICEN,1)
  364. CONS0.VAL(6,ICEN)=MPGAMN.VPOCHA(ICEN,1)
  365. ENDDO
  366. SEGDES MPRN
  367. SEGDES MPRETN
  368. SEGDES MPGN
  369. SEGDES MPGAMN
  370. C
  371. SEGINI, CONS1=CONS0
  372. C
  373. C**** Storage of l-1 values of F at the l Jacobi iteration
  374. C
  375. IPT1=MELTFA
  376. SEGACT IPT1
  377. NSOUPO= IPT1.LISOUS(/1)
  378. NSOUPO=MAX(NSOUPO,1)
  379. JG=NSOUPO
  380. SEGINI MLESUP
  381. IFACT=0
  382. IF(NSOUPO .EQ. 1)THEN
  383. MLESUP.LECT(1)=IPT1
  384. NFACE=IPT1.NUM(/1)
  385. NELT=IPT1.NUM(/2)
  386. IFACT=IFACT+(NFACE*NELT)
  387. ELSE
  388. DO ISOUP = 1, NSOUPO, 1
  389. IPT2=IPT1.LISOUS(ISOUP)
  390. MLESUP.LECT(ISOUP)=IPT2
  391. SEGACT IPT2
  392. NFACE=IPT2.NUM(/1)
  393. NELT=IPT2.NUM(/2)
  394. IFACT=IFACT+(NFACE*NELT)
  395. ENDDO
  396. ENDIF
  397. C
  398. C***********************************************************************
  399. C**** FIRST JACOBI ITERATION *******************************************
  400. C***********************************************************************
  401. C
  402. C
  403. C**** LOOP on ELTFA to compute:
  404. C
  405. C AN0_i = 1/(2 \Omega_i) \sum_j \sigma_i,j S_j
  406. C SIGMA0_j
  407. C The cut-off in each cell
  408. C The reciprocal of the dual time step USDTI
  409. C
  410. ICEN=0
  411. IFACG=0
  412. DO ISOUP=1,NSOUPO,1
  413. IPT2=MLESUP.LECT(ISOUP)
  414. NFACE=IPT2.NUM(/1)
  415. NELT=IPT2.NUM(/2)
  416. DO IELT=1,NELT,1
  417. ICEN=ICEN+1
  418. DO IFAC=1,NFACE,1
  419. C IFACG increase with IFAC,IELT,ISOUP
  420. IFACG = IFACG+1
  421. NGCF = IPT2.NUM(IFAC,IELT)
  422. NLCF = MLENT2.LECT(NGCF)
  423. C
  424. C************* NLCF = numero local du centre de facel
  425. C NGCF = numero global du centre de facel
  426. C NGCEG = numero global du centre ELT "gauche"
  427. C NLCEG = numero local du centre ELT "gauche"
  428. C NGCED = numero global du centre ELT "droite"
  429. C NLCED = numero local du centre ELT "droite"
  430. C
  431. NGCEG = MELEFL.NUM(1,NLCF)
  432. NGCED = MELEFL.NUM(3,NLCF)
  433. NLCEG = MLENT1.LECT(NGCEG)
  434. NLCED = MLENT1.LECT(NGCED)
  435. C
  436. CNX = MPNORM.VPOCHA(NLCF,7)
  437. CNY = MPNORM.VPOCHA(NLCF,8)
  438. CNZ = MPNORM.VPOCHA(NLCF,9)
  439. CTX = MPNORM.VPOCHA(NLCF,1)
  440. CTY = MPNORM.VPOCHA(NLCF,2)
  441. CTZ = MPNORM.VPOCHA(NLCF,3)
  442. CVX = MPNORM.VPOCHA(NLCF,4)
  443. CVY = MPNORM.VPOCHA(NLCF,5)
  444. CVZ = MPNORM.VPOCHA(NLCF,6)
  445. C
  446. C************* Left or right?
  447. C
  448. IF(NLCEG .EQ. ICEN)THEN
  449. LOGINV=.FALSE.
  450. ELSEIF(NLCED .EQ. ICEN)THEN
  451. LOGINV=.TRUE.
  452. ELSE
  453. WRITE(IOIMP,*) 'Anomalie dans kfm13.eso'
  454. WRITE(IOIMP,*) 'Probleme de SPG'
  455. C 21 2
  456. C Données incompatibles
  457. CALL ERREUR(21)
  458. GOTO 9999
  459. ENDIF
  460. C************* If right, inversion
  461. IF(LOGINV) THEN
  462. NLCED=NLCEG
  463. NLCEG=ICEN
  464. CNX = -1.0D0*CNX
  465. CNY = -1.0D0*CNY
  466. CNZ = -1.0D0*CNZ
  467. CTX = -1.0D0*CTX
  468. CTY = -1.0D0*CTY
  469. CTZ = -1.0D0*CTZ
  470. CVX = -1.0D0*CVX
  471. CVY = -1.0D0*CVY
  472. CVZ = -1.0D0*CVZ
  473. ENDIF
  474. C
  475. SURF = MPOVSU.VPOCHA(NLCF,1)
  476. C
  477. C************* Recuperation des Etats "gauche" et "droite"
  478. C
  479. C
  480. ROG = CONS0.VAL(1,NLCEG)
  481. RUXG = CONS0.VAL(2,NLCEG)
  482. RUYG = CONS0.VAL(3,NLCEG)
  483. RUZG = CONS0.VAL(4,NLCEG)
  484. RETG = CONS0.VAL(5,NLCEG)
  485. GAMG = CONS0.VAL(6,NLCEG)
  486. C
  487. UXG = RUXG / ROG
  488. UYG = RUYG / ROG
  489. UZG = RUZG / ROG
  490. C
  491. UNG = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ)
  492. UTG = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ)
  493. UVG = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ)
  494. C
  495. RECG = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
  496. RECG = RECG / ROG
  497. PG = (GAMG - 1.0D0) * (RETG - RECG)
  498. C
  499. VOLG = MPOVOL.VPOCHA(NLCEG,1)
  500. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  501. C
  502. ROD = CONS0.VAL(1,NLCED)
  503. RUXD = CONS0.VAL(2,NLCED)
  504. RUYD = CONS0.VAL(3,NLCED)
  505. RUZD = CONS0.VAL(4,NLCED)
  506. RETD = CONS0.VAL(5,NLCED)
  507. GAMD = CONS0.VAL(6,NLCED)
  508. C
  509. UXD = RUXD / ROD
  510. UYD = RUYD / ROD
  511. UZD = RUZD / ROD
  512. C
  513. UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
  514. UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
  515. UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
  516. C
  517. RECD = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2)
  518. RECD = RECD / ROD
  519. PD = (GAMD - 1.0D0) * (RETD - RECD)
  520. C
  521. IF(NLCEG .EQ. NLCED)THEN
  522. C Murs ou condition limite
  523. NLLIM=MLELIM.LECT(NGCF)
  524. IF(NLLIM .EQ.0)THEN
  525. C Mur
  526. UND = -1.0D0 * UNG
  527. UTD = UTG
  528. UVD = UVG
  529. UXD = UND * CNX + UTD * CTX + UVD * CVX
  530. UYD = UND * CNY + UTD * CTY + UVD * CVY
  531. UZD = UND * CNZ + UTD * CTZ + UVD * CVZ
  532. C
  533. RUXD = ROD * UXD
  534. RUYD = ROD * UYD
  535. RUZD = ROD * UZD
  536. C
  537. ELSE
  538. ROD = MPLIM.VPOCHA(NLLIM,1)
  539. UXD = MPLIM.VPOCHA(NLLIM,2)
  540. UYD = MPLIM.VPOCHA(NLLIM,3)
  541. UZD = MPLIM.VPOCHA(NLLIM,4)
  542. PD = MPLIM.VPOCHA(NLLIM,5)
  543. C
  544. RUXD = ROD*UXD
  545. RUYD = ROD*UYD
  546. RUZD = ROD * UZD
  547. C
  548. GAMD = GAMG
  549. C
  550. UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
  551. UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
  552. UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
  553. C
  554. RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
  555. & +(0.5D0*ROD*(UXD**2+UYD**2+UZD**2))
  556. ENDIF
  557. ENDIF
  558. C
  559. C************* We check that density and pressure are positive
  560. C
  561. IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  562. & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  563. WRITE(IOIMP,*) 'FMM'
  564. WRITE(IOIMP,*) 'TEST'
  565. C 22 0
  566. C**************** Opération malvenue. Résultat douteux
  567. C
  568. CALL ERREUR(22)
  569. IPROBL = 1
  570. GOTO 9998
  571. ENDIF
  572. C
  573. C************* We compute the interfacial state
  574. C
  575. ROF = 0.5D0*(ROG+ROD)
  576. PF = 0.5D0*(PG+PD)
  577. UNF = 0.5D0*(UNG+UND)
  578. GAMF = 0.5D0*(GAMG+GAMD)
  579. C
  580. C************* We compute SIGMA and the corrected cut-off speed
  581. C
  582. UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1))
  583. UCOF=MAX(UCOF,((UNG**2+UTG**2+UVG**2)**0.5D0))
  584. UCOF=MAX(UCOF,((UND**2+UTD**2+UVD**2)**0.5D0))
  585. IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1))
  586. & MPOCO1.VPOCHA(NLCEG,1)= UCOF
  587. C
  588. QN = ABS(UNF)
  589. CF = (GAMF * PF / ROF)**0.5D0
  590. IF(UCOF .GE. CF)THEN
  591. PHI2 = 1.0D0
  592. ELSEIF(QN .GT. CF)THEN
  593. PHI2 = 1.0D0
  594. ELSEIF(QN .LT. UCOF)THEN
  595. PHI2 = UCOF / CF
  596. ELSE
  597. PHI2 = QN / CF
  598. ENDIF
  599. PHI2 = PHI2 * PHI2
  600. C
  601. SIGMAF = (1.0D0 - PHI2) * QN
  602. SIGMAF = SIGMAF * SIGMAF
  603. SIGMAF = SIGMAF + (4.0D0 * PHI2 * CF * CF)
  604. SIGMAF = SIGMAF**0.5D0
  605. SIGMAF = SIGMAF + ((1.0D0 + PHI2) * QN)
  606. SIGMAF = 0.5D0 * SIGMAF
  607. C
  608. SIGMA0.VALVEC(NLCF) = SIGMAF
  609. C
  610. C************* We compute AN0
  611. C
  612. AN0.VALVEC(NLCEG) = AN0.VALVEC(NLCEG)
  613. & + (0.5D0 * SIGMAF * SURF / VOLG)
  614. C
  615. C************* We compute the dual time step
  616. C
  617. UNSDT = USDTI.VALVEC(NLCEG)
  618. UNSDTF = SIGMAF / DIAMG
  619. IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG) = UNSDTF
  620. C
  621. C************* We compute SIGMAV
  622. C
  623. C COEF=FG
  624. COEF = (MCOORD.XCOOR(((NGCEG-1)*4)+1)
  625. & -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX
  626. COEF = COEF
  627. & + (MCOORD.XCOOR(((NGCEG-1)*4)+2)
  628. & -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY
  629. COEF = COEF
  630. & + (MCOORD.XCOOR(((NGCEG-1)*4)+3)
  631. & -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ
  632. C
  633. C COEF1=FD
  634. COEF1 = (MCOORD.XCOOR(((NGCED-1)*4)+1)
  635. & -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX
  636. COEF1 = COEF1
  637. & + (MCOORD.XCOOR(((NGCED-1)*4)+2)
  638. & -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY
  639. COEF1 = COEF1
  640. & + (MCOORD.XCOOR(((NGCED-1)*4)+3)
  641. & -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ
  642. C
  643. C COEF=|FG|+|FD|
  644. COEF = ABS(COEF)+ABS(COEF1)
  645. SIGMAF = 0.5D0 * (MPCOEV.VPOCHA(NLCEG,1)
  646. & +MPCOEV.VPOCHA(NLCED,1)) / COEF
  647. C
  648. C
  649. SIGMAV.VALVEC(NLCF) = SIGMAF
  650. C
  651. C************* We compute the contribution of SIGMAV to
  652. C the diagonal part of the matrix to invert
  653. C
  654. C In the given report, DNO is a_i^0, i.e.
  655. C DN0 = AN0 (previously defined) +
  656. C + (\frac{1}{{\Omega_i}}
  657. C \sum_j \sigma_{V,i,j}^{0} S_j)
  658. C + 1 / \Delta t^n
  659. C
  660. COEF = SURF / VOLG
  661. DN0.VALVEC(NLCEG) = DN0.VALVEC(NLCEG) + COEF * SIGMAF
  662. C
  663. C
  664. C************* We compute the flux
  665. C FLU.VALVEC(1:4,IFAC) flux in IFAC face
  666. C 1 -> mass
  667. C 2 -> momentum (x)
  668. C 3 -> momentum (y)
  669. C 4 -> energy
  670. C
  671. RHTD = RETD + PD
  672. FR = (ROD * UND)
  673. FRUX = (ROD * UND * UXD) + (PD * CNX)
  674. FRUY = (ROD * UND * UYD) + (PD * CNY)
  675. FRUZ = (ROD * UND * UZD) + (PD * CNZ)
  676. FRET = (UND * RHTD)
  677. C
  678. COEF = 0.5D0 * SURF / VOLG
  679. C
  680. BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
  681. & + FR * COEF
  682. BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
  683. & + FRUX * COEF
  684. BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
  685. & + FRUY * COEF
  686. BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
  687. & + FRUZ * COEF
  688. BN0.VAL(5,NLCEG) = BN0.VAL(5,NLCEG)
  689. & + FRET * COEF
  690. C
  691. C************* Viscous contribution
  692. C
  693. COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLG
  694. C
  695. BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
  696. & - (COEF * ROD)
  697. BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
  698. & - (COEF * RUXD)
  699. BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
  700. & - (COEF * RUYD)
  701. BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
  702. & - (COEF * RUZD)
  703. BN0.VAL(5,NLCEG) = BN0.VAL(5,NLCEG)
  704. & - (COEF * RETD)
  705. C
  706. C************* On calcule CN0
  707. C
  708. C CN0 in the referred report is given by
  709. C \sum_j
  710. C \left\{
  711. C \frac{S_j}{2 \Omega_i} ~
  712. C \left(
  713. C \sigma_{i,j}^{0}
  714. C \mathbf{U}_{j,R}^{n+1,m}
  715. C \right)
  716. C \right\}
  717. C
  718. C
  719. COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
  720. C
  721. CN0.VAL(1,NLCEG) = CN0.VAL(1,NLCEG)
  722. & + (COEF * ROD)
  723. CN0.VAL(2,NLCEG) = CN0.VAL(2,NLCEG)
  724. & + (COEF * RUXD)
  725. CN0.VAL(3,NLCEG) = CN0.VAL(3,NLCEG)
  726. & + (COEF * RUYD)
  727. CN0.VAL(4,NLCEG) = CN0.VAL(4,NLCEG)
  728. & + (COEF * RUZD)
  729. CN0.VAL(5,NLCEG) = CN0.VAL(5,NLCEG)
  730. & + (COEF * RETD)
  731.  
  732. ENDDO
  733. C
  734. C********** End of loop on the faces of the element
  735. C
  736. C N.B.: ICEN = NLCEG
  737. C
  738. C
  739. C Computation of the low-Mach preconditioning Matrix
  740. C
  741. CG = (GAMG * PG / ROG)**0.5D0
  742. UCOG = MPOCO1.VPOCHA(ICEN,1)
  743. QG = 2.0D0 * RECG / 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. Q2.VAL(1,ICEN) = (0.5D0 * QG * QG)
  759. Q2.VAL(2,ICEN) = -1.0D0 * RUXG / ROG
  760. Q2.VAL(3,ICEN) = -1.0D0 * RUYG / ROG
  761. Q2.VAL(4,ICEN) = -1.0D0 * RUZG / ROG
  762. Q2.VAL(5,ICEN) = 1.0D0
  763. C
  764. COEF1 = (GAMG - 1.0D0) / (CG * CG)
  765. COEF = (1.0D0 / COEF1) + (0.5D0 * QG * QG)
  766. C COEF=HT
  767. Q1.VAL(1,ICEN) = COEF1
  768. Q1.VAL(2,ICEN) = COEF1 * RUXG / ROG
  769. Q1.VAL(3,ICEN) = COEF1 * RUYG / ROG
  770. Q1.VAL(4,ICEN) = COEF1 * RUZG / ROG
  771. Q1.VAL(5,ICEN) = COEF1 * COEF
  772. C
  773. C
  774. C Computation of the diagonal coefficients
  775. C
  776. AN0.VALVEC(ICEN) = AN0.VALVEC(ICEN)
  777. & + (USDTI.VALVEC(ICEN) * 2.0D0 / DCFL)
  778. C
  779. DN0.VALVEC(ICEN) = DN0.VALVEC(ICEN) + (1.0D0 / DTPS)
  780. c & + (USDTI.VALVEC(ICEN) * 2.0D0 / DCFL)
  781. C
  782. ENDDO
  783. ENDDO
  784. C
  785. C***********************************************************************
  786. C**** END FIRST JACOBI ITERATION ***************************************
  787. C***********************************************************************
  788. C
  789. C***********************************************************************
  790. C**** JACOBI LOOP ******************************************************
  791. C***********************************************************************
  792. C
  793. C
  794. LOGJAC=.TRUE.
  795. C
  796. DO IBLJAC=1,NJAC,1
  797. C
  798. C
  799. C Si ICOEF=2 (i.e. si on a "LJACOF" (voir kfm1.eso))
  800. C alors relaxation SGS avec uniquement des balayages "forward"
  801. C
  802. C Si ICOEF=3 (i.e. si on a "LJACOB")
  803. C alors relaxation SGS avec uniquement des balayages "backward"
  804. C
  805. C Sinon (i.e. si on a "LJACOFB")
  806. C alors relaxation SGS avec des balayages "forward" et "backward"
  807. C selon la parité de la sous-itération en cours (IBLJAC)...
  808. C
  809. IF(ICOEF .EQ. 2)THEN
  810. LOGJAC=.TRUE.
  811. ELSEIF(ICOEF .EQ. 3)THEN
  812. LOGJAC=.FALSE.
  813. ELSE
  814. LOGINV=LOGJAC
  815. LOGJAC=(IBLJAC / 2 * 2) .EQ. IBLJAC
  816. IF(LOGJAC)THEN
  817. LOGJAC= .NOT. LOGINV
  818. ELSE
  819. LOGJAC= LOGINV
  820. ENDIF
  821. ENDIF
  822. C
  823. IF(LOGJAC)THEN
  824. ICEN=0
  825. IFACG=0
  826. ID=1
  827. ISOUP1=1
  828. ISOUP2=NSOUPO
  829. ELSE
  830. ICEN=NCEN+1
  831. IFACG=IFACT+1
  832. ID=-1
  833. ISOUP1=NSOUPO
  834. ISOUP2=1
  835. ENDIF
  836. C
  837. DO ISOUP=ISOUP1,ISOUP2,ID
  838. IPT2=MLESUP.LECT(ISOUP)
  839. NFACE=IPT2.NUM(/1)
  840. NELT=IPT2.NUM(/2)
  841. C
  842. IF(LOGJAC)THEN
  843. IELT1=1
  844. IELT2=NELT
  845. IFAC1=1
  846. IFAC2=NFACE
  847. ELSE
  848. IELT2=1
  849. IELT1=NELT
  850. IFAC2=1
  851. IFAC1=NFACE
  852. ENDIF
  853. C
  854. DO IELT=IELT1,IELT2,ID
  855. ICEN=ICEN+ID
  856. DO IFAC=IFAC1,IFAC2,ID
  857. C IFACG increase with IFAC,IELT,ISOUP
  858. IFACG = IFACG+ID
  859. NGCF = IPT2.NUM(IFAC,IELT)
  860. NLCF = MLENT2.LECT(NGCF)
  861. NGCEG = MELEFL.NUM(1,NLCF)
  862. NGCED = MELEFL.NUM(3,NLCF)
  863. NLCEG = MLENT1.LECT(NGCEG)
  864. NLCED = MLENT1.LECT(NGCED)
  865. C
  866. CNX = MPNORM.VPOCHA(NLCF,7)
  867. CNY = MPNORM.VPOCHA(NLCF,8)
  868. CNZ = MPNORM.VPOCHA(NLCF,9)
  869. CTX = MPNORM.VPOCHA(NLCF,1)
  870. CTY = MPNORM.VPOCHA(NLCF,2)
  871. CTZ = MPNORM.VPOCHA(NLCF,3)
  872. CVX = MPNORM.VPOCHA(NLCF,4)
  873. CVY = MPNORM.VPOCHA(NLCF,5)
  874. CVZ = MPNORM.VPOCHA(NLCF,6)
  875. C
  876. C**************** Left or right?
  877. C
  878. IF(NLCEG .EQ. ICEN)THEN
  879. LOGINV=.FALSE.
  880. ELSEIF(NLCED .EQ. ICEN)THEN
  881. LOGINV=.TRUE.
  882. ELSE
  883. WRITE(IOIMP,*) 'Anomalie dans kfm14.eso'
  884. WRITE(IOIMP,*) 'Problème de SPG'
  885. C 21 2
  886. C Données incompatibles
  887. CALL ERREUR(21)
  888. GOTO 9999
  889. ENDIF
  890. C**************** If right, inversion
  891. IF(LOGINV) THEN
  892. NLCED=NLCEG
  893. NLCEG=ICEN
  894. CNX = -1.0D0*CNX
  895. CNY = -1.0D0*CNY
  896. CNZ = -1.0D0*CNZ
  897. CTX = -1.0D0*CTX
  898. CTY = -1.0D0*CTY
  899. CTZ = -1.0D0*CTZ
  900. CVX = -1.0D0*CVX
  901. CVY = -1.0D0*CVY
  902. CVZ = -1.0D0*CVZ
  903. ENDIF
  904. C
  905. SURF = MPOVSU.VPOCHA(NLCF,1)
  906. SIGMAF = SIGMA0.VALVEC(NLCF)
  907. C
  908. C**************** Recuperation des Etats "gauche" et "droite"
  909. C
  910. ROG = CONS1.VAL(1,NLCEG)
  911. RUXG = CONS1.VAL(2,NLCEG)
  912. RUYG = CONS1.VAL(3,NLCEG)
  913. RUZG = CONS1.VAL(4,NLCEG)
  914. RETG = CONS1.VAL(5,NLCEG)
  915. GAMG = CONS1.VAL(6,NLCEG)
  916. C
  917. UXG = RUXG / ROG
  918. UYG = RUYG / ROG
  919. UZG = RUZG / ROG
  920. C
  921. UNG = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ)
  922. UTG = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ)
  923. UVG = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ)
  924. C
  925. RECG = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
  926. RECG = RECG / ROG
  927. PG = (GAMG - 1.0D0) * (RETG - RECG)
  928. C
  929. VOLG = MPOVOL.VPOCHA(NLCEG,1)
  930. C
  931. ROD = CONS1.VAL(1,NLCED)
  932. RUXD = CONS1.VAL(2,NLCED)
  933. RUYD = CONS1.VAL(3,NLCED)
  934. RUZD = CONS1.VAL(4,NLCED)
  935. RETD = CONS1.VAL(5,NLCED)
  936. GAMD = CONS1.VAL(6,NLCED)
  937. C
  938. UXD = RUXD / ROD
  939. UYD = RUYD / ROD
  940. UZD = RUZD / ROD
  941. C
  942. UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
  943. UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
  944. UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
  945. C
  946. RECD = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2)
  947. RECD = RECD / ROD
  948. PD = (GAMD - 1.0D0) * (RETD - RECD)
  949. C
  950. IF(NLCEG .EQ. NLCED)THEN
  951. C Murs ou condition limite
  952. NLLIM=MLELIM.LECT(NGCF)
  953. IF(NLLIM .EQ.0)THEN
  954. C Mur
  955. UND = -1.0D0 * UNG
  956. UTD = UTG
  957. UVD = UVG
  958. UXD = UND * CNX + UTD * CTX + UVD * CVX
  959. UYD = UND * CNY + UTD * CTY + UVD * CVY
  960. UZD = UND * CNZ + UTD * CTZ + UVD * CVZ
  961. C
  962. RUXD = ROD*UXD
  963. RUYD = ROD*UYD
  964. RUZD = ROD * UZD
  965. C
  966. ELSE
  967. ROD = MPLIM.VPOCHA(NLLIM,1)
  968. UXD = MPLIM.VPOCHA(NLLIM,2)
  969. UYD = MPLIM.VPOCHA(NLLIM,3)
  970. UZD = MPLIM.VPOCHA(NLLIM,4)
  971. PD = MPLIM.VPOCHA(NLLIM,5)
  972. C
  973. RUXD = ROD*UXD
  974. RUYD = ROD*UYD
  975. RUZD = ROD * UZD
  976. C
  977. GAMD = GAMG
  978. C
  979. UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
  980. UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
  981. UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
  982. RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
  983. & +(0.5D0*ROD*(UXD**2+UYD**2+UZD**2))
  984. ENDIF
  985. ENDIF
  986. C
  987. C**************** We check that density and pressure are positive
  988. C
  989. IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
  990. & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
  991. WRITE(IOIMP,*) 'FMM'
  992. WRITE(IOIMP,*) 'ITERATION JAC',IBLJAC
  993. WRITE(IOIMP,*) 'ITERATION ICEN',ICEN
  994. WRITE(IOIMP,*) ROG,PG
  995. WRITE(IOIMP,*) ROD,PD
  996. C
  997. C 22 0
  998. C******************* Opération malvenue. Résultat douteux
  999. C
  1000. CALL ERREUR(22)
  1001. IPROBL = 1
  1002. GOTO 9998
  1003. ENDIF
  1004. C
  1005. C
  1006. C**************** We compute the flux
  1007. C FLU.VALVEC(1:4,IFAC) flux in IFAC face
  1008. C 1 -> mass
  1009. C 2 -> momentum (x)
  1010. C 3 -> momentum (y)
  1011. C 4 -> energy
  1012. C
  1013. RHTD = RETD + PD
  1014. FR = (ROD * UND)
  1015. FRUX = (ROD * UND * UXD) + (PD * CNX)
  1016. FRUY = (ROD * UND * UYD) + (PD * CNY)
  1017. FRUZ = (ROD * UND * UZD) + (PD * CNZ)
  1018. FRET = (UND * RHTD)
  1019. C
  1020. COEF = 0.5D0 * SURF / VOLG
  1021. C
  1022. BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  1023. & + FR * COEF
  1024. BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  1025. & + FRUX * COEF
  1026. BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  1027. & + FRUY * COEF
  1028. BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  1029. & + FRUZ * COEF
  1030. BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG)
  1031. & + FRET * COEF
  1032. C
  1033. C**************** Viscous contribution
  1034. C
  1035. COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLG
  1036. C
  1037. BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
  1038. & - (COEF * ROD)
  1039. BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
  1040. & - (COEF * RUXD)
  1041. BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
  1042. & - (COEF * RUYD)
  1043. BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
  1044. & - (COEF * RUZD)
  1045. BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG)
  1046. & - (COEF * RETD)
  1047. C
  1048. C**************** On calcule DIS1
  1049. C
  1050. C DIS1 in the referred report is given by
  1051. C \sum_j
  1052. C \left\{
  1053. C \frac{S_j}{2 \Omega_i} ~
  1054. C \left(
  1055. C \sigma_{i,j}^{0}
  1056. C \mathbf{U}_{j,R}^{n+1,m,l}
  1057. C \right)
  1058. C \right\}
  1059. C
  1060. C
  1061. COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
  1062. C
  1063. CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG)
  1064. & + (COEF * ROD)
  1065. CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG)
  1066. & + (COEF * RUXD)
  1067. CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG)
  1068. & + (COEF * RUYD)
  1069. CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG)
  1070. & + (COEF * RUZD)
  1071. CN1.VAL(5,NLCEG) = CN1.VAL(5,NLCEG)
  1072. & + (COEF * RETD)
  1073.  
  1074. ENDDO
  1075. C
  1076. C************* End of the loop on the faces of the element
  1077. C
  1078. DO ICOMP=1,5,1
  1079. D1.VALVEC(ICOMP) = 0.0D0
  1080. ENDDO
  1081. C
  1082. COEF=0.0D0
  1083. DO I1=1,5,1
  1084. COEF = COEF
  1085. & + (Q2.VAL(I1,ICEN) * (CN1.VAL(I1,ICEN)
  1086. & -CN0.VAL(I1,ICEN)))
  1087. ENDDO
  1088. COEF = COEF
  1089. & * (1.0D0 - PHI2V.VALVEC(ICEN))
  1090. & / PHI2V.VALVEC(ICEN)
  1091. C
  1092. DO I1=1,5,1
  1093. D1.VALVEC(I1) = (CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN))
  1094. & + (COEF * Q1.VAL(I1,ICEN))
  1095. ENDDO
  1096. C
  1097. C************* At this moment D1 = \Gamma \cdot (\Delta CN1)
  1098. C
  1099. DO ICOMP=1,5,1
  1100. C
  1101. C************* First increment
  1102. C
  1103. D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP)
  1104. & + RES.VAL(ICOMP,ICEN)
  1105. & + (BN0.VAL(ICOMP,ICEN)
  1106. & -BN1.VAL(ICOMP,ICEN)))
  1107. & / (AN0.VALVEC(ICEN)
  1108. & +DN0.VALVEC(ICEN))
  1109. ENDDO
  1110. C
  1111. C************* At this moment, in the referred report,
  1112. C D1 is (\delta \mathbf{U}')^{l+1}_i in section A3
  1113. C
  1114. C************* Second increment
  1115. C
  1116. COEF=0.0D0
  1117. DO I1=1,5,1
  1118. COEF = COEF + (Q2.VAL(I1,ICEN) * D1.VALVEC(I1))
  1119. ENDDO
  1120. C
  1121. C************* COEF is the scalar product between Q2 and the first
  1122. C increment.
  1123. C Let us multiply it by
  1124. C b_i^0 / (a_i^0 + b_i^0)
  1125. C
  1126. COEF = COEF
  1127. & * AN0.VALVEC(ICEN)
  1128. & * (PHI2V.VALVEC(ICEN) - 1.0D0)
  1129. COEF = COEF
  1130. & / (AN0.VALVEC(ICEN)
  1131. & + (DN0.VALVEC(ICEN) * PHI2V.VALVEC(ICEN)))
  1132. C
  1133. DO I1=1,5,1
  1134. D1.VALVEC(I1) = D1.VALVEC(I1)
  1135. & + COEF * Q1.VAL(I1,ICEN)
  1136. ENDDO
  1137. C
  1138. C************* D1 = \delta \mathbf{U}
  1139. C
  1140. DO ICOMP=1,5,1
  1141. DU.VAL(ICOMP,ICEN) = D1.VALVEC(ICOMP)
  1142. C
  1143. CONS1.VAL(ICOMP,ICEN) = CONS0.VAL(ICOMP,ICEN)
  1144. & + D1.VALVEC(ICOMP)
  1145. C
  1146. BN1.VAL(ICOMP,ICEN) = 0.0D0
  1147. CN1.VAL(ICOMP,ICEN) = 0.0D0
  1148. ENDDO
  1149. C
  1150. C********** End of the loop on the elements
  1151. ENDDO
  1152. C
  1153. C******* End of the loop on the Domains
  1154. ENDDO
  1155. C
  1156. C**** End of the Jacobi loop
  1157. ENDDO
  1158. C
  1159. 9998 CONTINUE
  1160. DO ICEN=1,NCEN,1
  1161. DO ICOMP=1,5,1
  1162. MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
  1163. ENDDO
  1164. ENDDO
  1165. C
  1166. SEGSUP DU
  1167. SEGSUP RES
  1168. SEGSUP CONS0
  1169. SEGSUP CONS1
  1170. SEGSUP Q1
  1171. SEGSUP Q2
  1172. SEGSUP SIGMA0
  1173. SEGSUP AN0
  1174. SEGSUP AN0
  1175. SEGSUP BN0
  1176. SEGSUP BN1
  1177. SEGSUP CN0
  1178. SEGSUP CN1
  1179. SEGSUP DN0
  1180. SEGSUP PHI2V
  1181. SEGSUP USDTI
  1182. C
  1183. SEGDES MPOVSU
  1184. SEGDES MPOVDI
  1185. SEGDES MPOVOL
  1186. SEGDES MPNORM
  1187. SEGDES MPOCO
  1188. SEGDES MPCOEV
  1189. C
  1190. SEGSUP MPOCO1
  1191. C
  1192. SEGDES MPODU
  1193. C
  1194. SEGDES MELEFL
  1195. SEGSUP MLENT1
  1196. SEGSUP MLENT2
  1197. C
  1198. SEGSUP MLELIM
  1199. IF(MPLIM .GT. 0) SEGDES MPLIM
  1200. C
  1201. DO ISOUP=1,NSOUPO,1
  1202. IPT2=MLESUP.LECT(ISOUP)
  1203. SEGDES IPT2
  1204. ENDDO
  1205. IPT2 = MELTFA
  1206. IF(ISOUP .GT. 1) SEGDES IPT2
  1207. SEGSUP MLESUP
  1208. C
  1209. 9999 CONTINUE
  1210. RETURN
  1211. END
  1212. C
  1213.  
  1214.  
  1215.  
  1216.  

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