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

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