Télécharger kgfm12.eso

Retour à la liste

Numérotation des lignes :

  1. C KGFM12 SOURCE BECC 11/05/26 21:15:48 6981
  2. SUBROUTINE KGFM12(
  3. & LOGCON, NESP,
  4. & INDMET,
  5. & MLRMGA, MLRMPI, MLRPGA, MLRPPI,
  6. & PMIN,
  7. & ICHPHI, ICHALP,
  8. & IPHI1,IROF1,IVITF1,IPF1,
  9. & IFRMAF, IFRALF,
  10. & ICHPSU,ICHPDI,ICHPVO,
  11. & MELEMC,MELEMF,MELEFE,MELLIM,
  12. & ICHRES,
  13. & DT,
  14. & LOGNC,LOGAN,MESERR)
  15. C************************************************************************
  16. C
  17. C PROJET : CASTEM 2000
  18. C
  19. C NOM : KGFM12
  20. C
  21. C DESCRIPTION : Voir KONV15
  22. C
  23. C Cas deux dimensions, GFMP
  24. C
  25. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  26. C
  27. C AUTEUR : A. BECCANTINI, DEN/DM2S/SFME/LTMF
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES
  32. C
  33. C 1) Parameters concerning some gas properties
  34. C
  35. C GAM1, PINF1, GAM2, PINF2 : gamma and pinf for the "stiffened gases".
  36. C
  37. C 2) Numerical scheme
  38. C
  39. C INDMET : 1 GODUNOV (exact Riemann solver)
  40. C
  41. C 3) Pointers for the CHPOINTs/MCHAML des variables
  42. C
  43. C ICHPHI : chpoint CENTRE contenant (phi)
  44. C
  45. C IPHI1 : MCHAML sur "FACEL" contenant phi
  46. C ("gauche" et "droite");
  47. C
  48. C IROF1 : MCHAML sur "FACEL" contenant la masse volumique 1
  49. C ("gauche" et "droite");
  50. C
  51. C IVITF1 : MCHAML sur "FACEL" contenant la vitesse 1 dans le repaire
  52. C local (n,t) et les cosinus directeurs des repaire local;
  53. C
  54. C IPF1 : MCHAML sur "FACEL" contenant la pression 1;
  55. C
  56. C
  57. C 4) Others
  58. C
  59. C ICHPSU : CHPOINT "FACE" contenant la surface des faces
  60. C
  61. C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
  62. C de chaque element
  63. C
  64. C ICHPVO : CHPOINT "CENTRE" contenant le volume
  65. C de chaque element
  66. C
  67. C 5) Pointeurs de MELEME de la table DOMAINE
  68. C
  69. C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
  70. C
  71. C MELEMF : MELEME 'FACE' du SPG des FACES
  72. C
  73. C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
  74. C
  75. C MELLIM : MAILLAGE where the flux (or residu) will not be found
  76. C
  77. C SORTIES (il faudrait dire E/S)
  78. C
  79. C ICHRES : pointeurs de CHPOINTs "FACE" des flux aux interfaces:
  80. C
  81. C DT : pas de temps pour le respect de la CFL-like condition
  82. C DT < DIAMEL /2 /max(Lambda_i)
  83. C En maillage regulier cette condition garantie la
  84. C non-interaction des ondes
  85. C
  86. C LOGNC : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
  87. C dans pour la solution du probleme Riemann n'a pas bien
  88. C marchéee
  89. C
  90. C LOGAN : (LOGICAL): si .TRUE. une anomalie à été detectée
  91. C
  92. C MESERR : pour l'ecriture des messages d'erreurs
  93. C
  94. C************************************************************************
  95. C
  96. C HISTORIQUE (Anomalies et modifications éventuelles)
  97. C
  98. C 04/12/2010 - Created
  99. * 25/05/2011 - Evolution in CAST3M
  100. C
  101. C************************************************************************
  102. C
  103. C
  104. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  105. C Y \in (0,1)
  106. C Si non il faut le faire!!!
  107. C
  108. C************************************************************************
  109. C
  110. C
  111. IMPLICIT INTEGER(I-N)
  112. INTEGER NESP, NESP1, IESP, INDMET, NPAR
  113. & , ICHPHI, ICHALP
  114. & , IPHI1, IROF1, IVITF1, IPF1
  115. & , ICHPSU, ICHPDI, ICHPVO, MELEMC, MELEMF, MELEFE, MELLIM
  116. & , IGEOMC, IGEOMF
  117. & , ICHRES
  118. & , NFAC
  119. & , NLCF, NGCEG, NGCED, NLCEG, NLCED
  120. & , NGCF, NLCF1, SPG1, SPG2, NLFL
  121. & , NFRA
  122. REAL*8 GAMMA, PINF, ALP, ALPI, RNUM, DEN, PMIN
  123. REAL*8 DT, UNSDT, CELLT, CELLTM
  124. & , PHIG1C, PHIG1F, ROG1, UNG1, UTG1, PG1
  125. & , PHID1C, PHID1F, ROD1, UND1, UTD1, PD1
  126. & , SURF,CNX, CNY, CTX , CTY
  127. & , CELL, DIAMG, DIAMD, DIAM, VOLG, VOLD
  128. & , FLUCEL, FLUCED
  129. & , PHIPRO, PINFG, GAMG, PINFD, GAMD
  130. REAL*8 RHO_B, P_B, U_B, D_B
  131. & ,RHO_A, P_A, U_A, D_A
  132. & ,RHO_B1, P_B1, U_B1, D_B1
  133. & ,RHO_A1, P_A1, U_A1, D_A1
  134. & ,ROG1G, PG1G, UNG1G
  135. & ,ROD1G, PD1G, UND1G
  136. & ,UINT
  137. C
  138. LOGICAL LOGNC, LOGAN, LOGVAC, LOGDEB, LOGCON
  139. CHARACTER*(40) MESERR
  140. CHARACTER*(8) TYPE
  141. PARAMETER (LOGDEB = .FALSE.)
  142. C
  143. C**** LES INCLUDES
  144. C
  145. -INC CCOPTIO
  146. -INC SMCOORD
  147. -INC SMELEME
  148. -INC SMCHAML
  149. C Normal and tangent
  150. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
  151. & MELT1X.MELVAL, MELT1Y.MELVAL
  152. POINTEUR MELUN1.MELVAL, MELUT1.MELVAL
  153. POINTEUR MELPH1.MELVAL, MELRO1.MELVAL, MELP1.MELVAL
  154. POINTEUR MCHAMY.MCHAML, MCHAMA.MCHAML
  155. -INC SMCHPOI
  156. POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL, MPOVVO.MPOVAL
  157. & , MPORES.MPOVAL
  158. & , MPPHI1.MPOVAL, MPALPH.MPOVAL
  159. -INC SMLMOTS
  160. -INC SMLENTI
  161. POINTEUR MLELIM.MLENTI
  162. INTEGER JG
  163. -INC SMLREEL
  164. POINTEUR MLRMGA.MLREEL, MLRPGA.MLREEL,
  165. & MLRMPI.MLREEL, MLRPPI.MLREEL
  166. C
  167. C**** Segment for the subroutine for the flux
  168. C
  169. POINTEUR IFLU.MLREEL, WPARG.MLREEL, WPARD.MLREEL
  170. C
  171. C**** Pointer for Y
  172. SEGMENT FRAMAS
  173. REAL*8 YET(NFRA)
  174. ENDSEGMENT
  175. POINTEUR FRAMAG.FRAMAS, FRAMAD.FRAMAS
  176. & , ALPHAG.FRAMAS, ALPHAD.FRAMAS
  177. & , ALPHCG.FRAMAS, ALPHCD.FRAMAS
  178. C
  179. LOGAN = .FALSE.
  180. C
  181. C**** Parameter (phi + fractions massiques)
  182. C
  183. NPAR = 1 + NESP
  184. JG = NPAR
  185. SEGINI WPARG
  186. SEGINI WPARD
  187. C
  188. SEGACT MLRMGA
  189. SEGACT MLRPGA
  190. SEGACT MLRMPI
  191. SEGACT MLRPPI
  192. C
  193. C**** KRIPAD pour la correspondance global/local des conditions limits
  194. C
  195. CALL KRIPAD(MELLIM,MLELIM)
  196. C SEGINI MLELIM
  197. C
  198. C**** Initialisation des MCHAMLs
  199. C
  200. C PHI
  201. C
  202. MCHEL1 = IPHI1
  203. SEGACT MCHEL1
  204. MCHAM1 = MCHEL1.ICHAML(1)
  205. SEGACT MCHAM1
  206. MELPH1 = MCHAM1.IELVAL(1)
  207. SEGDES MCHEL1
  208. SEGDES MCHAM1
  209. C
  210. C Masse volumique
  211. C
  212. MCHEL1 = IROF1
  213. SEGACT MCHEL1
  214. MCHAM1 = MCHEL1.ICHAML(1)
  215. SEGACT MCHAM1
  216. MELRO1 = MCHAM1.IELVAL(1)
  217. SEGDES MCHEL1
  218. SEGDES MCHAM1
  219. C
  220. C Pression
  221. C
  222. MCHEL1 = IPF1
  223. SEGACT MCHEL1
  224. MCHAM1 = MCHEL1.ICHAML(1)
  225. SEGACT MCHAM1
  226. MELP1 = MCHAM1.IELVAL(1)
  227. SEGDES MCHEL1
  228. SEGDES MCHAM1
  229. C
  230. C Vitesse et cosinus directeurs du repere (n,t)
  231. C
  232. MCHEL1 = IVITF1
  233. SEGACT MCHEL1
  234. C
  235. C La vitesse a comme SPG MELEFE
  236. C Le cosinus directeurs ont comme SPG MELEMF
  237. C
  238. C MCHAM1 -> Cosinus directeurs
  239. C MCHAM2 -> Vitesse
  240. C
  241. SPG1 = MCHEL1.IMACHE(1)
  242. SPG2 = MCHEL1.IMACHE(2)
  243. IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
  244. MCHAM1 = MCHEL1.ICHAML(1)
  245. MCHAM2 = MCHEL1.ICHAML(2)
  246. ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
  247. MCHAM1 = MCHEL1.ICHAML(2)
  248. MCHAM2 = MCHEL1.ICHAML(1)
  249. ELSE
  250. LOGAN = .TRUE.
  251. CALL ERREUR(251)
  252. GOTO 9999
  253. ENDIF
  254. SEGACT MCHAM1
  255. MELVNX = MCHAM1.IELVAL(1)
  256. MELVNY = MCHAM1.IELVAL(2)
  257. MELT1X = MCHAM1.IELVAL(3)
  258. MELT1Y = MCHAM1.IELVAL(4)
  259. SEGDES MCHAM1
  260. SEGACT MCHAM2
  261. MELUN1 = MCHAM2.IELVAL(1)
  262. MELUT1 = MCHAM2.IELVAL(2)
  263. SEGDES MCHAM2
  264. SEGDES MCHEL1
  265. C
  266. IF (NESP .GE. 1) THEN
  267. C
  268. C Fractions massiques au centre (non-conservative approach)
  269. C
  270. CALL LICHT(ICHALP, MPALPH, TYPE, IGEOMC)
  271. C SEGACT MPALPH
  272.  
  273. NFRA = NESP
  274. SEGINI FRAMAG
  275. SEGINI FRAMAD
  276. SEGINI ALPHAG
  277. SEGINI ALPHAD
  278. SEGINI ALPHCG
  279. SEGINI ALPHCD
  280. C
  281. MCHEL1 = IFRMAF
  282. SEGACT MCHEL1
  283. MCHAMY = MCHEL1.ICHAML(1)
  284. SEGACT MCHAMY
  285. SEGDES MCHEL1
  286. NESP1 = MCHAMY.IELVAL(/1)
  287. IF (NESP1 .NE. NESP) THEN
  288. WRITE(IOIMP,*) 'Mass fractions, dimension = ???'
  289. CALL ERREUR(21)
  290. GOTO 9999
  291. ENDIF
  292. C
  293. DO IESP = 1, NESP, 1
  294. MELVA1 = MCHAMY.IELVAL(IESP)
  295. SEGACT MELVA1
  296. ENDDO
  297. C
  298. MCHEL1 = IFRALF
  299. SEGACT MCHEL1
  300. MCHAMA = MCHEL1.ICHAML(1)
  301. SEGACT MCHAMA
  302. SEGDES MCHEL1
  303. NESP1 = MCHAMA.IELVAL(/1)
  304. IF (NESP1 .NE. NESP) THEN
  305. WRITE(IOIMP,*) 'Volume fractions dimension = ???'
  306. CALL ERREUR(21)
  307. GOTO 9999
  308. ENDIF
  309. C
  310. DO IESP = 1, NESP, 1
  311. MELVA1 = MCHAMA.IELVAL(IESP)
  312. SEGACT MELVA1
  313. * write(*,*) melva1
  314. ENDDO
  315. ENDIF
  316. C
  317. C write(*,*) 'Son in kgdm12'
  318. C write(*,*) 'Ho letto gli MCHAM'
  319. C Fluxes
  320. C
  321. JG = IDIM + 3 + NESP
  322. SEGINI IFLU
  323. C
  324. C**** Initialisation des MELEMEs
  325. C
  326. C 'CENTRE', 'FACEL'
  327. C
  328. IPT2 = MELEFE
  329. SEGACT IPT2
  330. NFAC = IPT2.NUM(/2)
  331. C
  332. C**** KRIPAD pour la correspondance global/local de centre
  333. C
  334. CALL KRIPAD(MELEMC,MLENT1)
  335. C
  336. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  337. C
  338. C Si i est le numero global d'un noeud de ICEN,
  339. C MLENT1.LECT(i) contient sa position, i.e.
  340. C
  341. C I = numero global du noeud centre
  342. C MLENT1.LECT(i) = numero local du noeud centre
  343. C
  344. C MLENT1 déjà activé, i.e.
  345. C
  346. C SEGINI MLENT1
  347. C
  348. C
  349. C**** KRIPAD pour la correspondance global/local de 'FACE'
  350. C
  351. CALL KRIPAD(MELEMF,MLENT2)
  352. C SEGINI MLENT2
  353. C
  354. C**** CHPOINTs de la table DOMAINE
  355. C
  356. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  357. CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  358. CALL LICHT(ICHPVO,MPOVVO,TYPE,IGEOMC)
  359. C
  360. C**** LICHT active les MPOVALs en *MOD
  361. C
  362. C i.e.
  363. C
  364. C SEGACT MPOVSU*MOD
  365. C SEGACT MPOVDI*MOD
  366. C SEGACT MPOVVO*MOD
  367. C
  368. C
  369. C**** Les FLUX aux face
  370. C
  371. CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
  372. C SEGACT MPORES*MOD
  373. C
  374. C**** Activation des MCHAMLs
  375. C
  376. SEGACT MELPH1
  377. SEGACT MELRO1
  378. SEGACT MELP1
  379. SEGACT MELUN1
  380. SEGACT MELUT1
  381. SEGACT MELVNX
  382. SEGACT MELVNY
  383. SEGACT MELT1X
  384. SEGACT MELT1Y
  385. C
  386. C**** Lecture des MPOVALs ICHPHI, ICHPA2
  387. C
  388. CALL LICHT(ICHPHI, MPPHI1, TYPE, IGEOMC)
  389. C SEGACT MPPHI1
  390. C
  391. C**** Initialisation de 1/DT
  392. C
  393. UNSDT = 0.0D0
  394. C
  395. C**** BOUCLE SUR FACEL pour le calcul du FLUX/RESIDU
  396. C
  397. DO NLCF = 1, NFAC
  398. C
  399. C******* NLCF = numero local du centre de facel
  400. C NGCF = numero global du centre de facel
  401. C NLCF1 = numero local du centre de face
  402. C NGCEG = numero global du centre ELT "gauche"
  403. C NLCEG = numero local du centre ELT "gauche"
  404. C NGCED = numero global du centre ELT "droite"
  405. C NLCED = numero local du centre ELT "droite"
  406. C
  407. NGCEG = IPT2.NUM(1,NLCF)
  408. NGCED = IPT2.NUM(3,NLCF)
  409. NGCF = IPT2.NUM(2,NLCF)
  410. NLCF1 = MLENT2.LECT(NGCF)
  411. NLCEG = MLENT1.LECT(NGCEG)
  412. NLCED = MLENT1.LECT(NGCED)
  413. NLFL = MLELIM.LECT(NGCF)
  414. C
  415. VOLG = MPOVVO.VPOCHA(NLCEG,1)
  416. VOLD = MPOVVO.VPOCHA(NLCED,1)
  417. C
  418. C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
  419. C
  420. IF(NLCF .NE. NLCF1)THEN
  421. MESERR = 'Il ne faut pas jouer avec la console. '
  422. LOGAN = .TRUE.
  423. CALL ERREUR(5)
  424. GOTO 9999
  425. ENDIF
  426. IF(NLFL .EQ. 0)THEN
  427. CELLTM = 0.0D0
  428. C
  429. C********** Recuperation des Etats "gauche" et "droite"
  430. C
  431. PHIG1C = MPPHI1.VPOCHA(NLCEG,1)
  432. PHIG1F = MELPH1.VELCHE(1,NLCF)
  433. ROG1 = MELRO1.VELCHE(1,NLCF)
  434. UNG1 = MELUN1.VELCHE(1,NLCF)
  435. UTG1 = MELUT1.VELCHE(1,NLCF)
  436. PG1 = MELP1.VELCHE(1,NLCF)
  437. DO IESP = 1, NESP, 1
  438. MELVA1 = MCHAMY.IELVAL(IESP)
  439. FRAMAG.YET(IESP) = MELVA1.VELCHE(1,NLCF)
  440. ENDDO
  441. DO IESP = 1, NESP, 1
  442. MELVA1 = MCHAMA.IELVAL(IESP)
  443. ALPHAG.YET(IESP) = MELVA1.VELCHE(1,NLCF)
  444. C IF (NGCF .EQ. 685) THEN
  445. C write(*,*) 'IFRALF =', IFRALF
  446. C write(*,*) 'MELVA1 =', MELVA1
  447. C write(*,*) 'NLCF =', NLCF
  448. C write(*,*) 'IESP =', IESP
  449. C write(*,*) 'ALPHAG.YET(IESP) =', MELVA1.VELCHE(1,NLCF)
  450. C ENDIF
  451. ENDDO
  452. DO IESP = 1, NESP, 1
  453. ALPHCG . YET(IESP) = MPALPH.VPOCHA(NLCEG,IESP)
  454. ENDDO
  455. C
  456. WPARG.PROG(1) = PHIG1F
  457. C Fractions massiques traité en forme conservative...
  458. DO IESP = 1, NESP, 1
  459. WPARG.PROG (IESP + 1) = FRAMAG.YET(IESP)
  460. ENDDO
  461. C
  462. PHID1C = MPPHI1.VPOCHA(NLCED,1)
  463. PHID1F = MELPH1.VELCHE(3,NLCF)
  464. ROD1 = MELRO1.VELCHE(3,NLCF)
  465. UND1 = MELUN1.VELCHE(3,NLCF)
  466. UTD1 = MELUT1.VELCHE(3,NLCF)
  467. PD1 = MELP1.VELCHE(3,NLCF)
  468. DO IESP = 1, NESP, 1
  469. MELVA1 = MCHAMY.IELVAL(IESP)
  470. FRAMAD.YET(IESP) = MELVA1.VELCHE(3,NLCF)
  471. ENDDO
  472. DO IESP = 1, NESP, 1
  473. MELVA1 = MCHAMA.IELVAL(IESP)
  474. ALPHAD.YET(IESP) = MELVA1.VELCHE(3,NLCF)
  475. ENDDO
  476. DO IESP = 1, NESP, 1
  477. ALPHCD . YET(IESP) = MPALPH.VPOCHA(NLCED,IESP)
  478. ENDDO
  479. WPARD.PROG(1) = PHID1F
  480. DO IESP = 1, NESP, 1
  481. WPARD.PROG (IESP + 1) = FRAMAD.YET(IESP)
  482. ENDDO
  483. C
  484. CNX = MELVNX.VELCHE(1,NLCF)
  485. CNY = MELVNY.VELCHE(1,NLCF)
  486. CTX = MELT1X.VELCHE(1,NLCF)
  487. CTY = MELT1Y.VELCHE(1,NLCF)
  488. C
  489. SURF = MPOVSU.VPOCHA(NLCF,1)
  490. C
  491. PHIPRO = PHIG1C * PHID1C
  492. C
  493. C IF (NGCF .EQ. 685) THEN
  494. IF (LOGDEB) THEN
  495. WRITE(*,*)
  496. WRITE(*,*) 'NGCF =', NGCF
  497. WRITE(*,*) 'Left '
  498. WRITE(*,*) 'Phic, phi'
  499. WRITE(*,*) PHIG1C, PHIG1F
  500. WRITE(*,*) 'rho, un, ut, pg'
  501. WRITE(*,*) ROG1, UNG1, UTG1, PG1
  502. WRITE(*,*) 'Right'
  503. WRITE(*,*) 'Phic, phi'
  504. WRITE(*,*) PHID1C, PHID1F
  505. WRITE(*,*) 'rho, un, ut, pg'
  506. WRITE(*,*) ROD1, UND1, UTD1, PD1
  507. ENDIF
  508. C
  509. C********** According to PHIG1C we compute GAMG, PINFG
  510. C
  511. IF (PHIG1C .LT. 0.0D0) THEN
  512. MLREE1 = MLRMGA
  513. MLREE2 = MLRMPI
  514. ELSE
  515. MLREE1 = MLRPGA
  516. MLREE2 = MLRPPI
  517. ENDIF
  518. C Left
  519. ALP = 1.0D0
  520. DEN = 0.0D0
  521. RNUM = 0.0D0
  522. DO IESP = 1, NESP, 1
  523. GAMMA = MLREE1.PROG(IESP)
  524. PINF = MLREE2.PROG(IESP)
  525. ALPI = ALPHAG.YET(IESP)
  526. C IF (NGCF .EQ. 685) THEN
  527. C write(*,*) 'IESP =', IESP
  528. C write(*,*) 'ALPI =', ALPI
  529. C ENDIF
  530. C
  531. ALP = ALP - ALPI
  532. DEN = DEN + (ALPI / (GAMMA - 1.0D0))
  533. RNUM = RNUM +
  534. & ((ALPI * GAMMA * PINF) / (GAMMA - 1.0D0))
  535. C IF (NGCF .EQ. 685) THEN
  536. C write(*,*) 'IESP =', IESP
  537. C write(*,*) 'GAMMA, PINF, ALPI, ALP, DEN, RNUM'
  538. C write(*,*) GAMMA, PINF, ALPI, ALP, DEN, RNUM
  539. C ENDIF
  540. ENDDO
  541. GAMMA = MLREE1.PROG(NESP + 1)
  542. PINF = MLREE2.PROG(NESP + 1)
  543. DEN = DEN + (ALP / (GAMMA - 1.0D0))
  544. RNUM = RNUM + ((ALP * GAMMA * PINF) / (GAMMA - 1.0D0))
  545. C IF (NGCF .EQ. 685) THEN
  546. C write(*,*) 'IESP =', (NESP + 1)
  547. C write(*,*) 'GAMMA, PINF, ALPI, ALP, DEN, RNUM'
  548. C write(*,*) GAMMA, PINF, ALPI, ALP, DEN, RNUM
  549. C ENDIF
  550. C
  551. PINF = RNUM / DEN
  552. GAMMA = 1.0D0 / DEN
  553. GAMMA = GAMMA + 1.0D0
  554. PINF = PINF / GAMMA
  555. C
  556. PINFG = PINF
  557. GAMG = GAMMA
  558. C
  559. C********** According to PHID1C we compute GAMD, PINFD
  560. C
  561. IF (PHID1C .LT. 0.0D0) THEN
  562. MLREE1 = MLRMGA
  563. MLREE2 = MLRMPI
  564. ELSE
  565. MLREE1 = MLRPGA
  566. MLREE2 = MLRPPI
  567. ENDIF
  568. ALP = 1.0D0
  569. DEN = 0.0D0
  570. RNUM = 0.0D0
  571. DO IESP = 1, NESP, 1
  572. GAMMA = MLREE1.PROG(IESP)
  573. PINF = MLREE2.PROG(IESP)
  574. ALPI = ALPHAD.YET(IESP)
  575. ALP = ALP - ALPI
  576. DEN = DEN + (ALPI / (GAMMA - 1.0D0))
  577. RNUM = RNUM +
  578. & ((ALPI * GAMMA * PINF) / (GAMMA - 1.0D0))
  579. ENDDO
  580. GAMMA = MLREE1.PROG(NESP + 1)
  581. PINF = MLREE2.PROG(NESP + 1)
  582. DEN = DEN + (ALP / (GAMMA - 1.0D0))
  583. RNUM = RNUM + ((ALP * GAMMA * PINF) / (GAMMA - 1.0D0))
  584. C
  585. PINF = RNUM / DEN
  586. GAMMA = 1.0D0 / DEN
  587. GAMMA = GAMMA + 1.0D0
  588. PINF = PINF / GAMMA
  589. C
  590. PINFD = PINF
  591. GAMD = GAMMA
  592. C
  593. C A priori, there is no vacuum...
  594. C
  595. LOGVAC = .FALSE.
  596. C
  597. C IF (NGCF .EQ. 685) THEN
  598. C WRITE(*,*)
  599. C & PINFG, GAMG,
  600. C & PINFD, GAMD
  601. C ENDIF
  602. C
  603. IF (INDMET .EQ. 1) THEN
  604. C
  605. C Exact Riemann solver
  606. C Computation of the states on the left and on the
  607. C right of the contact discontinuity
  608. C
  609. CALL FSTIFF(IDIM, NPAR,
  610. & PINFG, GAMG,
  611. & PINFD, GAMD,
  612. & PMIN,
  613. & ROG1, PG1, UNG1, UTG1, 0.0D0, WPARG.PROG,
  614. & ROD1, PD1, UND1, UTD1, 0.0D0, WPARD.PROG,
  615. & RHO_B, P_B, U_B, D_B,
  616. & RHO_A, P_A, U_A, D_A,
  617. & 0.0D0,
  618. & UINT,
  619. & IFLU.PROG, CELLT,
  620. & LOGAN, LOGNC, LOGVAC)
  621. IF (LOGAN) THEN
  622. WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  623. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  624. C
  625. WRITE(*,*)
  626. WRITE(*,*) 'NGCF =', NGCF
  627. WRITE(*,*) 'NLCF =', NLCF
  628. WRITE(*,*) 'Left '
  629. WRITE(*,*) 'Phic, phi'
  630. WRITE(*,*) PHIG1C, PHIG1F
  631. WRITE(*,*) 'Pinf, gamma'
  632. WRITE(*,*) PINFG, GAMG
  633. WRITE(*,*) 'rho, un, ut, pg'
  634. WRITE(*,*) ROG1, UNG1, UTG1, PG1
  635. WRITE(*,*) 'Right'
  636. WRITE(*,*) 'Phic, phi'
  637. WRITE(*,*) PHID1C, PHID1F
  638. WRITE(*,*) 'Pinf, gamma'
  639. WRITE(*,*) PINFD, GAMD
  640. WRITE(*,*) 'rho, un, ut, pg'
  641. WRITE(*,*) ROD1, UND1, UTD1, PD1
  642. C
  643. CALL ERREUR(5)
  644. GOTO 9999
  645. ENDIF
  646. ENDIF
  647. CELLTM = MAX(CELLT, CELLTM)
  648. C
  649. IF ((.NOT. LOGVAC) .AND. (PHIPRO .GT. 0.0D0)) THEN
  650. C
  651. C************* Ecriture du residu
  652. C
  653. C FLUX(2) = RO Un RO Un
  654. C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
  655. C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
  656. C FLUX(5) = RO Un Et RO Un Et
  657. C Mass fractions... (5 + IESP)...
  658. C Volume fractions in non conservative form...
  659. C
  660. FLUCEL = IFLU.PROG(1) * SURF
  661. C write(*,*) 'flux', iflu.prog(1)
  662. C write(*,*) 'nlceg, nlced', nlceg, nlced
  663. C write(*,*) 'surf', surf
  664. MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
  665. & (FLUCEL / VOLG)
  666. IF (NLCED .NE. NLCEG)
  667. & MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
  668. & (FLUCEL / VOLD)
  669. FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
  670. MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
  671. & (FLUCEL / VOLG)
  672. IF (NLCED .NE. NLCEG)
  673. & MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
  674. & (FLUCEL / VOLD)
  675. FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
  676. MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
  677. & (FLUCEL / VOLG)
  678. IF (NLCED .NE. NLCEG)
  679. & MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
  680. & (FLUCEL / VOLD)
  681. FLUCEL = IFLU.PROG(4) * SURF
  682. MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
  683. & (FLUCEL / VOLG)
  684. IF (NLCED .NE. NLCEG)
  685. & MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
  686. & (FLUCEL / VOLD)
  687. C
  688. C************* rho Y in conservative form
  689. C Components 5 + 1 -> 5 + NESP
  690. C
  691. DO IESP = 1, NESP, 1
  692.  
  693. FLUCEL = IFLU.PROG(5 + IESP) * SURF
  694. MPORES.VPOCHA(NLCEG,5 + IESP) =
  695. & MPORES.VPOCHA(NLCEG,5 + IESP) -
  696. & (FLUCEL / VOLG)
  697. IF (NLCED .NE. NLCEG)
  698. & MPORES.VPOCHA(NLCED,5 + IESP) =
  699. & MPORES.VPOCHA(NLCED,5 + IESP) +
  700. & (FLUCEL / VOLD)
  701. ENDDO
  702. C************* Alpha in non conservative form
  703. C Components 5 + NESP -> 5 + (2 * NESP)
  704. IF (NESP .GE. 1) THEN
  705. IF (UINT .GT. 0.0D0) THEN
  706. FRAMAS = ALPHAG
  707. ELSE
  708. FRAMAS = ALPHAD
  709. ENDIF
  710. DO IESP = 1, NESP, 1
  711. C UINT > 0
  712. C FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
  713. C UINT < 0
  714. C FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
  715. FLUCEL = UINT * (FRAMAS.YET(IESP) -
  716. & ALPHCG.YET(IESP)) * SURF
  717. C UINT > 0
  718. C FLUCED = UINT * (PHIG1F - PHID1C) * SURF
  719. C UINT < 0
  720. C FLUCED = UINT * (PHID1F - PHID1C) * SURF
  721. FLUCED = UINT * (FRAMAS.YET(IESP) -
  722. & ALPHCD.YET(IESP)) * SURF
  723. C
  724. MPORES.VPOCHA(NLCEG,5 + NESP + IESP) =
  725. & MPORES.VPOCHA(NLCEG,5 + NESP + IESP) -
  726. & (FLUCEL / VOLG)
  727. IF (NLCED .NE. NLCEG)
  728. & MPORES.VPOCHA(NLCED,5 + NESP + IESP) =
  729. & MPORES.VPOCHA(NLCED,5 + NESP + IESP) +
  730. & (FLUCED / VOLD)
  731. ENDDO
  732. ENDIF
  733. ENDIF
  734. C
  735. C Level set
  736. C
  737. IF (LOGCON) THEN
  738. C
  739. C************* Conservative, variable rho phi...
  740. C
  741. FLUCEL = IFLU.PROG(5) * SURF
  742. C write(*,*) 'flux', iflu.prog(1)
  743. C write(*,*) 'nlceg, nlced', nlceg, nlced
  744. C write(*,*) 'surf', surf
  745. MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
  746. & (FLUCEL / VOLG)
  747. IF (NLCED .NE. NLCEG)
  748. & MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
  749. & (FLUCEL / VOLD)
  750. ELSE
  751. C
  752. C************ Non conservative, variable phi
  753. C
  754. IF (UINT .GT. 0.0D0) THEN
  755. FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
  756. FLUCED = UINT * (PHIG1F - PHID1C) * SURF
  757. ELSE
  758. FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
  759. FLUCED = UINT * (PHID1F - PHID1C) * SURF
  760. ENDIF
  761. C
  762. MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
  763. & (FLUCEL / VOLG)
  764. IF (NLCED .NE. NLCEG)
  765. & MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
  766. & (FLUCED / VOLD)
  767. ENDIF
  768. C
  769. IF ((PHIPRO .LE. 0.0D0) .OR. LOGVAC) THEN
  770. * write(*,*) 'Son qui !!!'
  771. * stop
  772. C
  773. C Left flux for everything but not for the level-set
  774. C
  775. ROD1G = RHO_B
  776. PD1G = P_B
  777. UND1G = U_B
  778. IF (INDMET .EQ. 1) THEN
  779. C
  780. C Exact Riemann solver
  781. C
  782. CALL FSTIFF(IDIM, NPAR,
  783. & PINFG, GAMG,
  784. & PINFG, GAMG,
  785. & PMIN,
  786. & ROG1, PG1, UNG1, UTG1, 0.0D0, WPARG.PROG,
  787. & ROD1G, PD1G, UND1G, UTG1, 0.0D0, WPARG.PROG,
  788. & RHO_B1, P_B1, U_B1, D_B1,
  789. & RHO_A1, P_A1, U_A1, D_A1,
  790. & 0.0D0,
  791. & UINT,
  792. & IFLU.PROG, CELLT,
  793. & LOGAN, LOGNC, LOGVAC)
  794. C IF (LOGVAC) THEN
  795. C WRITE(IOIMP,*) 'Vacuum of vacuum ???'
  796. C WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  797. C WRITE(IOIMP,*) 'ANOMALY DETECTED'
  798. C CALL ERREUR(251)
  799. C GOTO 9999
  800. C ENDIF
  801. IF (LOGAN) THEN
  802. WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  803. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  804. C
  805. WRITE(*,*)
  806. WRITE(*,*) 'NGCF =', NGCF
  807. WRITE(*,*) 'Left '
  808. WRITE(*,*) 'Phic, phi'
  809. WRITE(*,*) PHIG1C, PHIG1F
  810. WRITE(*,*) 'rho, un, ut, pg'
  811. WRITE(*,*) ROG1, UNG1, UTG1, PG1
  812. WRITE(*,*) 'Right'
  813. WRITE(*,*) 'Phic, phi'
  814. WRITE(*,*) PHID1C, PHID1F
  815. WRITE(*,*) 'rho, un, ut, pg'
  816. WRITE(*,*) ROD1, UND1, UTD1, PD1
  817. C
  818. CALL ERREUR(5)
  819. GOTO 9999
  820. ENDIF
  821. ENDIF
  822. CELLTM = MAX(CELLT, CELLTM)
  823. C
  824. C********** Ecriture du residu... gauche only...
  825. C
  826. C FLUX(2) = RO Un RO Un
  827. C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
  828. C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
  829. C FLUX(5) = RO Un Et RO Un Et
  830. C
  831. FLUCEL = IFLU.PROG(1) * SURF
  832. C write(*,*) 'flux', iflu.prog(1)
  833. C write(*,*) 'nlceg, nlced', nlceg, nlced
  834. C write(*,*) 'surf', surf
  835. MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
  836. & (FLUCEL / VOLG)
  837. FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
  838. MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
  839. & (FLUCEL / VOLG)
  840. FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
  841. MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
  842. & (FLUCEL / VOLG)
  843. FLUCEL = IFLU.PROG(4) * SURF
  844. MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
  845. & (FLUCEL / VOLG)
  846. C
  847. C Level set already done...
  848. C Level set treated as the other variables (like in the
  849. C following comment... does not work)
  850. C
  851. C FLUCEL = IFLU.PROG(5) * SURF
  852. C MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
  853. C & (FLUCEL / VOLG)
  854. C
  855. C
  856. C************* rho Y in conservative form
  857. C Components 5 + 1 -> 5 + NESP
  858. C
  859. DO IESP = 1, NESP, 1
  860. FLUCEL = IFLU.PROG(5 + IESP) * SURF
  861. MPORES.VPOCHA(NLCEG,5 + IESP) =
  862. & MPORES.VPOCHA(NLCEG,5 + IESP) -
  863. & (FLUCEL / VOLG)
  864. ENDDO
  865. C************* Alpha in non conservative form
  866. C Components 5 + NESP -> 5 + (2 * NESP)
  867. IF (NESP .GE. 1) THEN
  868. IF (UINT .GT. 0.0D0) THEN
  869. FRAMAS = ALPHAG
  870. ELSE
  871. FRAMAS = ALPHAD
  872. ENDIF
  873. DO IESP = 1, NESP, 1
  874. C UINT > 0
  875. C FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
  876. C UINT < 0
  877. C FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
  878. FLUCEL = UINT * (FRAMAS.YET(IESP) -
  879. & ALPHCG.YET(IESP)) * SURF
  880. C
  881. MPORES.VPOCHA(NLCEG,5 + NESP + IESP) =
  882. & MPORES.VPOCHA(NLCEG,5 + NESP + IESP) -
  883. & (FLUCEL / VOLG)
  884. ENDDO
  885. ENDIF
  886. C
  887. C********** Calcul du pas du temps (CFL)
  888. C
  889. C********* a) etat a l'interface
  890. C
  891. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  892. DIAMD = MPOVDI.VPOCHA(NLCED,1)
  893. DIAM = MIN(DIAMG,DIAMD)
  894. CELL = CELLTM/DIAM
  895. IF(CELL .GT. UNSDT)THEN
  896. UNSDT = CELL
  897. ENDIF
  898. C
  899. C Right flux for everything but not for the level-set
  900. C
  901. ROG1G = RHO_A
  902. PG1G = P_A
  903. UNG1G = U_A
  904. C
  905. IF (NLCED .NE. NLCEG) THEN
  906. IF (INDMET .EQ. 1) THEN
  907. C
  908. C Exact Riemann solver
  909. C
  910. CALL FSTIFF(IDIM, NPAR,
  911. & PINFD, GAMD,
  912. & PINFD, GAMD,
  913. & PMIN,
  914. & ROG1G, PG1G, UNG1G, UTD1, 0.0D0, WPARD.PROG,
  915. & ROD1, PD1, UND1, UTD1, 0.0D0, WPARD.PROG,
  916. & RHO_B1, P_B1, U_B1, D_B1,
  917. & RHO_A1, P_A1, U_A1, D_A1,
  918. & 0.0D0,
  919. & UINT,
  920. & IFLU.PROG, CELLT,
  921. & LOGAN, LOGNC, LOGVAC)
  922. C
  923. C IF (LOGVAC) THEN
  924. C WRITE(IOIMP,*) 'Vacuum of vacuum ???'
  925. C WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  926. C WRITE(IOIMP,*) 'ANOMALY DETECTED'
  927. C CALL ERREUR(251)
  928. C GOTO 9999
  929. C ENDIF
  930. C
  931. IF (LOGAN) THEN
  932. WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  933. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  934. CALL ERREUR(5)
  935. GOTO 9999
  936. ENDIF
  937. C
  938. ENDIF
  939. CELLTM = MAX(CELLT, CELLTM)
  940.  
  941. C
  942. C************* Ecriture du flux a droite...
  943. C
  944. C FLUX(2) = RO Un RO Un
  945. C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
  946. C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
  947. C FLUX(5) = RO Un Et RO Un Et
  948. C
  949. FLUCEL = IFLU.PROG(1) * SURF
  950. C write(*,*) 'flux', iflu.prog(1)
  951. C write(*,*) 'nlceg, nlced', nlceg, nlced
  952. C write(*,*) 'surf', surf
  953. MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
  954. & (FLUCEL / VOLD)
  955. FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
  956.  
  957. MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
  958. & (FLUCEL / VOLD)
  959. FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
  960. MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
  961. & (FLUCEL / VOLD)
  962. FLUCEL = IFLU.PROG(4) * SURF
  963. MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
  964. & (FLUCEL / VOLD)
  965. C
  966. C Level set already done...
  967. C Level set treated as the other variables (like in the
  968. C following comment... does not work)
  969. C
  970. C FLUCEL = IFLU.PROG(5) * SURF
  971. C MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
  972. C & (FLUCEL / VOLD)
  973. C
  974. C
  975. C**************** rho Y in conservative form
  976. C Components 5 + 1 -> 5 + NESP
  977. C
  978. DO IESP = 1, NESP, 1
  979.  
  980. FLUCEL = IFLU.PROG(5 + IESP) * SURF
  981. MPORES.VPOCHA(NLCED,5 + IESP) =
  982. & MPORES.VPOCHA(NLCED,5 + IESP) +
  983. & (FLUCEL / VOLD)
  984. ENDDO
  985. C************* Alpha in non conservative form
  986. C Components 5 + NESP -> 5 + (2 * NESP)
  987. IF (NESP .GE. 1) THEN
  988. IF (UINT .GT. 0.0D0) THEN
  989. FRAMAS = ALPHAG
  990. ELSE
  991. FRAMAS = ALPHAD
  992. ENDIF
  993. DO IESP = 1, NESP, 1
  994. FLUCED = UINT * (FRAMAS.YET(IESP) -
  995. & ALPHCD.YET(IESP)) * SURF
  996. C
  997. MPORES.VPOCHA(NLCED,5 + NESP + IESP) =
  998. & MPORES.VPOCHA(NLCED,5 + NESP + IESP) +
  999. & (FLUCED / VOLD)
  1000. ENDDO
  1001. ENDIF
  1002. ENDIF
  1003. ENDIF
  1004. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  1005. DIAMD = MPOVDI.VPOCHA(NLCED,1)
  1006. DIAM = MIN(DIAMG,DIAMD)
  1007. CELL = CELLTM/DIAM
  1008. IF(CELL .GT. UNSDT)THEN
  1009. UNSDT = CELL
  1010. ENDIF
  1011. ENDIF
  1012. C
  1013. C**** Fin boucle sur FACEL
  1014. C
  1015. ENDDO
  1016. C
  1017. C**** Pas du temps (condition de non interaction en 1D)
  1018. C
  1019. DT = 0.5D0 / UNSDT
  1020. C
  1021. C**** Desactivation des segments
  1022. C
  1023. SEGSUP MLELIM
  1024. SEGSUP MLENT1
  1025. SEGSUP MLENT2
  1026. C
  1027. SEGDES MPOVSU
  1028. SEGDES MPOVDI
  1029. SEGDES MPOVVO
  1030. SEGDES MPORES
  1031. SEGDES IPT2
  1032. C
  1033. C Warning. MLRMFR = PGAS . 'MASSFRA'
  1034. C MLRCHE = PGAS . 'CHEMCOEF'
  1035. C
  1036. SEGSUP WPARG
  1037. SEGSUP WPARD
  1038. C
  1039. SEGSUP IFLU
  1040. C
  1041. SEGDES MPPHI1
  1042. SEGDES MELPH1
  1043. SEGDES MELRO1
  1044. SEGDES MELP1
  1045. SEGDES MELUN1
  1046. SEGDES MELUT1
  1047. SEGDES MELVNX
  1048. SEGDES MELVNY
  1049. SEGDES MELT1X
  1050. SEGDES MELT1Y
  1051. C
  1052. SEGDES MLRMGA
  1053. SEGDES MLRPGA
  1054. SEGDES MLRMPI
  1055. SEGDES MLRPPI
  1056. C
  1057. IF (NESP .GE. 1)THEN
  1058. C
  1059. SEGDES MPALPH
  1060. C
  1061. SEGSUP FRAMAG
  1062. SEGSUP FRAMAD
  1063. SEGSUP ALPHAG
  1064. SEGSUP ALPHAD
  1065. SEGSUP ALPHCG
  1066. SEGSUP ALPHCD
  1067. C
  1068. DO IESP = 1, NESP
  1069. MELVA1 = MCHAMY.IELVAL(IESP)
  1070. SEGDES MELVA1
  1071. ENDDO
  1072. SEGDES MCHAMY
  1073. C
  1074. DO IESP = 1, NESP
  1075. MELVA1 = MCHAMA.IELVAL(IESP)
  1076. SEGDES MELVA1
  1077. ENDDO
  1078. SEGDES MCHAMA
  1079. ENDIF
  1080. C
  1081. 9999 CONTINUE
  1082. RETURN
  1083. END
  1084. CC
  1085.  
  1086.  
  1087.  
  1088.  

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