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
194 : CALL KRIPAD(MELLIM,MLELIM)
195 : C SEGINI MLELIM
196 : C
197 : C**** Initialisation des MCHAMLs
198 : C
199 : C PHI
200 : C
201 : MCHEL1 = IPHI1
202 : SEGACT MCHEL1
203 : MCHAM1 = MCHEL1.ICHAML(1)
204 : SEGACT MCHAM1
205 : MELPH1 = MCHAM1.IELVAL(1)
206 : SEGDES MCHEL1
207 : SEGDES MCHAM1
208 : C
209 : C Masse volumique
210 : C
211 : MCHEL1 = IROF1
212 : SEGACT MCHEL1
213 : MCHAM1 = MCHEL1.ICHAML(1)
214 : SEGACT MCHAM1
215 : MELRO1 = MCHAM1.IELVAL(1)
216 : SEGDES MCHEL1
217 : SEGDES MCHAM1
218 : C
219 : C Pression
220 : C
221 : MCHEL1 = IPF1
222 : SEGACT MCHEL1
223 : MCHAM1 = MCHEL1.ICHAML(1)
224 : SEGACT MCHAM1
225 : MELP1 = MCHAM1.IELVAL(1)
226 : SEGDES MCHEL1
227 : SEGDES MCHAM1
228 : C
229 : C Vitesse et cosinus directeurs du repere (n,t)
230 : C
231 : MCHEL1 = IVITF1
232 : SEGACT MCHEL1
233 : C
234 : C La vitesse a comme SPG MELEFE
235 : C Le cosinus directeurs ont comme SPG MELEMF
236 : C
237 : C MCHAM1 -> Cosinus directeurs
238 : C MCHAM2 -> Vitesse
239 : C
240 : SPG1 = MCHEL1.IMACHE(1)
241 : SPG2 = MCHEL1.IMACHE(2)
242 : IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
243 : MCHAM1 = MCHEL1.ICHAML(1)
244 : MCHAM2 = MCHEL1.ICHAML(2)
245 : ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
246 : MCHAM1 = MCHEL1.ICHAML(2)
247 : MCHAM2 = MCHEL1.ICHAML(1)
248 : ELSE
249 : LOGAN = .TRUE.
249 : CALL ERREUR(251)
250 : GOTO 9999
251 : ENDIF
252 : SEGACT MCHAM1
253 : MELVNX = MCHAM1.IELVAL(1)
254 : MELVNY = MCHAM1.IELVAL(2)
255 : MELT1X = MCHAM1.IELVAL(3)
256 : MELT1Y = MCHAM1.IELVAL(4)
257 : SEGDES MCHAM1
258 : SEGACT MCHAM2
259 : MELUN1 = MCHAM2.IELVAL(1)
260 : MELUT1 = MCHAM2.IELVAL(2)
261 : SEGDES MCHAM2
262 : SEGDES MCHEL1
263 : C
264 : IF (NESP .GE. 1) THEN
265 : C
266 : C Fractions massiques au centre (non-conservative approach)
267 : C
267 : CALL LICHT(ICHALP, MPALPH, TYPE, IGEOMC)
268 : C SEGACT MPALPH
269 :
270 : NFRA = NESP
271 : SEGINI FRAMAG
272 : SEGINI FRAMAD
273 : SEGINI ALPHAG
274 : SEGINI ALPHAD
275 : SEGINI ALPHCG
276 : SEGINI ALPHCD
277 : C
278 : MCHEL1 = IFRMAF
279 : SEGACT MCHEL1
280 : MCHAMY = MCHEL1.ICHAML(1)
281 : SEGACT MCHAMY
282 : SEGDES MCHEL1
283 : NESP1 = MCHAMY.IELVAL(/1)
284 : IF (NESP1 .NE. NESP) THEN
285 : WRITE(IOIMP,*) 'Mass fractions, dimension = ???'
285 : CALL ERREUR(21)
286 : GOTO 9999
287 : ENDIF
288 : C
289 : DO IESP = 1, NESP, 1
290 : MELVA1 = MCHAMY.IELVAL(IESP)
291 : SEGACT MELVA1
292 : ENDDO
293 : C
294 : MCHEL1 = IFRALF
295 : SEGACT MCHEL1
296 : MCHAMA = MCHEL1.ICHAML(1)
297 : SEGACT MCHAMA
298 : SEGDES MCHEL1
299 : NESP1 = MCHAMA.IELVAL(/1)
300 : IF (NESP1 .NE. NESP) THEN
301 : WRITE(IOIMP,*) 'Volume fractions dimension = ???'
301 : CALL ERREUR(21)
302 : GOTO 9999
303 : ENDIF
304 : C
305 : DO IESP = 1, NESP, 1
306 : MELVA1 = MCHAMA.IELVAL(IESP)
307 : SEGACT MELVA1
308 : * write(*,*) melva1
309 : ENDDO
310 : ENDIF
311 : C
312 : C write(*,*) 'Son in kgdm12'
313 : C write(*,*) 'Ho letto gli MCHAM'
314 : C Fluxes
315 : C
316 : JG = IDIM + 3 + NESP
317 : SEGINI IFLU
318 : C
319 : C**** Initialisation des MELEMEs
320 : C
321 : C 'CENTRE', 'FACEL'
322 : C
323 : IPT2 = MELEFE
324 : SEGACT IPT2
325 : NFAC = IPT2.NUM(/2)
326 : C
327 : C**** KRIPAD pour la correspondance global/local de centre
328 : C
328 : CALL KRIPAD(MELEMC,MLENT1)
329 : C
330 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
331 : C
332 : C Si i est le numero global d'un noeud de ICEN,
333 : C MLENT1.LECT(i) contient sa position, i.e.
334 : C
335 : C I = numero global du noeud centre
336 : C MLENT1.LECT(i) = numero local du noeud centre
337 : C
338 : C MLENT1 déjà activé, i.e.
339 : C
340 : C SEGINI MLENT1
341 : C
342 : C
343 : C**** KRIPAD pour la correspondance global/local de 'FACE'
344 : C
344 : CALL KRIPAD(MELEMF,MLENT2)
345 : C SEGINI MLENT2
346 : C
347 : C**** CHPOINTs de la table DOMAINE
348 : C
348 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
348 : CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
348 : CALL LICHT(ICHPVO,MPOVVO,TYPE,IGEOMC)
349 : C
350 : C**** LICHT active les MPOVALs en *MOD
351 : C
352 : C i.e.
353 : C
354 : C SEGACT MPOVSU*MOD
355 : C SEGACT MPOVDI*MOD
356 : C SEGACT MPOVVO*MOD
357 : C
358 : C
359 : C**** Les FLUX aux face
360 : C
360 : CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
361 : C SEGACT MPORES*MOD
362 : C
363 : C**** Activation des MCHAMLs
364 : C
365 : SEGACT MELPH1
366 : SEGACT MELRO1
367 : SEGACT MELP1
368 : SEGACT MELUN1
369 : SEGACT MELUT1
370 : SEGACT MELVNX
371 : SEGACT MELVNY
372 : SEGACT MELT1X
373 : SEGACT MELT1Y
374 : C
375 : C**** Lecture des MPOVALs ICHPHI, ICHPA2
376 : C
376 : CALL LICHT(ICHPHI, MPPHI1, TYPE, IGEOMC)
377 : C SEGACT MPPHI1
378 : C
379 : C**** Initialisation de 1/DT
380 : C
381 : UNSDT = 0.0D0
382 : C
383 : C**** BOUCLE SUR FACEL pour le calcul du FLUX/RESIDU
384 : C
385 : DO NLCF = 1, NFAC
386 : C
387 : C******* NLCF = numero local du centre de facel
388 : C NGCF = numero global du centre de facel
389 : C NLCF1 = numero local du centre de face
390 : C NGCEG = numero global du centre ELT "gauche"
391 : C NLCEG = numero local du centre ELT "gauche"
392 : C NGCED = numero global du centre ELT "droite"
393 : C NLCED = numero local du centre ELT "droite"
394 : C
395 : NGCEG = IPT2.NUM(1,NLCF)
396 : NGCED = IPT2.NUM(3,NLCF)
397 : NGCF = IPT2.NUM(2,NLCF)
398 : NLCF1 = MLENT2.LECT(NGCF)
399 : NLCEG = MLENT1.LECT(NGCEG)
400 : NLCED = MLENT1.LECT(NGCED)
401 : NLFL = MLELIM.LECT(NGCF)
402 : C
403 : VOLG = MPOVVO.VPOCHA(NLCEG,1)
404 : VOLD = MPOVVO.VPOCHA(NLCED,1)
405 : C
406 : C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
407 : C
408 : IF(NLCF .NE. NLCF1)THEN
409 : MESERR = 'Il ne faut pas jouer avec la console. '
410 : LOGAN = .TRUE.
410 : CALL ERREUR(5)
411 : GOTO 9999
412 : ENDIF
413 : IF(NLFL .EQ. 0)THEN
414 : CELLTM = 0.0D0
415 : C
416 : C********** Recuperation des Etats "gauche" et "droite"
417 : C
418 : PHIG1C = MPPHI1.VPOCHA(NLCEG,1)
419 : PHIG1F = MELPH1.VELCHE(1,NLCF)
420 : ROG1 = MELRO1.VELCHE(1,NLCF)
421 : UNG1 = MELUN1.VELCHE(1,NLCF)
422 : UTG1 = MELUT1.VELCHE(1,NLCF)
423 : PG1 = MELP1.VELCHE(1,NLCF)
424 : DO IESP = 1, NESP, 1
425 : MELVA1 = MCHAMY.IELVAL(IESP)
426 : FRAMAG.YET(IESP) = MELVA1.VELCHE(1,NLCF)
427 : ENDDO
428 : DO IESP = 1, NESP, 1
429 : MELVA1 = MCHAMA.IELVAL(IESP)
430 : ALPHAG.YET(IESP) = MELVA1.VELCHE(1,NLCF)
431 : C IF (NGCF .EQ. 685) THEN
432 : C write(*,*) 'IFRALF =', IFRALF
433 : C write(*,*) 'MELVA1 =', MELVA1
434 : C write(*,*) 'NLCF =', NLCF
435 : C write(*,*) 'IESP =', IESP
436 : C write(*,*) 'ALPHAG.YET(IESP) =', MELVA1.VELCHE(1,NLCF)
437 : C ENDIF
438 : ENDDO
439 : DO IESP = 1, NESP, 1
440 : ALPHCG . YET(IESP) = MPALPH.VPOCHA(NLCEG,IESP)
441 : ENDDO
442 : C
443 : WPARG.PROG(1) = PHIG1F
444 : C Fractions massiques traité en forme conservative...
445 : DO IESP = 1, NESP, 1
446 : WPARG.PROG (IESP + 1) = FRAMAG.YET(IESP)
447 : ENDDO
448 : C
449 : PHID1C = MPPHI1.VPOCHA(NLCED,1)
450 : PHID1F = MELPH1.VELCHE(3,NLCF)
451 : ROD1 = MELRO1.VELCHE(3,NLCF)
452 : UND1 = MELUN1.VELCHE(3,NLCF)
453 : UTD1 = MELUT1.VELCHE(3,NLCF)
454 : PD1 = MELP1.VELCHE(3,NLCF)
455 : DO IESP = 1, NESP, 1
456 : MELVA1 = MCHAMY.IELVAL(IESP)
457 : FRAMAD.YET(IESP) = MELVA1.VELCHE(3,NLCF)
458 : ENDDO
459 : DO IESP = 1, NESP, 1
460 : MELVA1 = MCHAMA.IELVAL(IESP)
461 : ALPHAD.YET(IESP) = MELVA1.VELCHE(3,NLCF)
462 : ENDDO
463 : DO IESP = 1, NESP, 1
464 : ALPHCD . YET(IESP) = MPALPH.VPOCHA(NLCED,IESP)
465 : ENDDO
466 : WPARD.PROG(1) = PHID1F
467 : DO IESP = 1, NESP, 1
468 : WPARD.PROG (IESP + 1) = FRAMAD.YET(IESP)
469 : ENDDO
470 : C
471 : CNX = MELVNX.VELCHE(1,NLCF)
472 : CNY = MELVNY.VELCHE(1,NLCF)
473 : CTX = MELT1X.VELCHE(1,NLCF)
474 : CTY = MELT1Y.VELCHE(1,NLCF)
475 : C
476 : SURF = MPOVSU.VPOCHA(NLCF,1)
477 : C
478 : PHIPRO = PHIG1C * PHID1C
479 : C
480 : C IF (NGCF .EQ. 685) THEN
481 : IF (LOGDEB) THEN
482 : WRITE(*,*)
483 : WRITE(*,*) 'NGCF =', NGCF
484 : WRITE(*,*) 'Left '
485 : WRITE(*,*) 'Phic, phi'
486 : WRITE(*,*) PHIG1C, PHIG1F
487 : WRITE(*,*) 'rho, un, ut, pg'
488 : WRITE(*,*) ROG1, UNG1, UTG1, PG1
489 : WRITE(*,*) 'Right'
490 : WRITE(*,*) 'Phic, phi'
491 : WRITE(*,*) PHID1C, PHID1F
492 : WRITE(*,*) 'rho, un, ut, pg'
493 : WRITE(*,*) ROD1, UND1, UTD1, PD1
494 : ENDIF
495 : C
496 : C********** According to PHIG1C we compute GAMG, PINFG
497 : C
498 : IF (PHIG1C .LT. 0.0D0) THEN
499 : MLREE1 = MLRMGA
500 : MLREE2 = MLRMPI
501 : ELSE
502 : MLREE1 = MLRPGA
503 : MLREE2 = MLRPPI
504 : ENDIF
505 : C Left
506 : ALP = 1.0D0
507 : DEN = 0.0D0
508 : RNUM = 0.0D0
509 : DO IESP = 1, NESP, 1
510 : GAMMA = MLREE1.PROG(IESP)
511 : PINF = MLREE2.PROG(IESP)
512 : ALPI = ALPHAG.YET(IESP)
513 : C IF (NGCF .EQ. 685) THEN
514 : C write(*,*) 'IESP =', IESP
515 : C write(*,*) 'ALPI =', ALPI
516 : C ENDIF
517 : C
518 : ALP = ALP - ALPI
519 : DEN = DEN + (ALPI / (GAMMA - 1.0D0))
520 : RNUM = RNUM +
521 : & ((ALPI * GAMMA * PINF) / (GAMMA - 1.0D0))
522 : C IF (NGCF .EQ. 685) THEN
523 : C write(*,*) 'IESP =', IESP
524 : C write(*,*) 'GAMMA, PINF, ALPI, ALP, DEN, RNUM'
525 : C write(*,*) GAMMA, PINF, ALPI, ALP, DEN, RNUM
526 : C ENDIF
527 : ENDDO
528 : GAMMA = MLREE1.PROG(NESP + 1)
529 : PINF = MLREE2.PROG(NESP + 1)
530 : DEN = DEN + (ALP / (GAMMA - 1.0D0))
531 : RNUM = RNUM + ((ALP * GAMMA * PINF) / (GAMMA - 1.0D0))
532 : C IF (NGCF .EQ. 685) THEN
533 : C write(*,*) 'IESP =', (NESP + 1)
534 : C write(*,*) 'GAMMA, PINF, ALPI, ALP, DEN, RNUM'
535 : C write(*,*) GAMMA, PINF, ALPI, ALP, DEN, RNUM
536 : C ENDIF
537 : C
538 : PINF = RNUM / DEN
539 : GAMMA = 1.0D0 / DEN
540 : GAMMA = GAMMA + 1.0D0
541 : PINF = PINF / GAMMA
542 : C
543 : PINFG = PINF
544 : GAMG = GAMMA
545 : C
546 : C********** According to PHID1C we compute GAMD, PINFD
547 : C
548 : IF (PHID1C .LT. 0.0D0) THEN
549 : MLREE1 = MLRMGA
550 : MLREE2 = MLRMPI
551 : ELSE
552 : MLREE1 = MLRPGA
553 : MLREE2 = MLRPPI
554 : ENDIF
555 : ALP = 1.0D0
556 : DEN = 0.0D0
557 : RNUM = 0.0D0
558 : DO IESP = 1, NESP, 1
559 : GAMMA = MLREE1.PROG(IESP)
560 : PINF = MLREE2.PROG(IESP)
561 : ALPI = ALPHAD.YET(IESP)
562 : ALP = ALP - ALPI
563 : DEN = DEN + (ALPI / (GAMMA - 1.0D0))
564 : RNUM = RNUM +
565 : & ((ALPI * GAMMA * PINF) / (GAMMA - 1.0D0))
566 : ENDDO
567 : GAMMA = MLREE1.PROG(NESP + 1)
568 : PINF = MLREE2.PROG(NESP + 1)
569 : DEN = DEN + (ALP / (GAMMA - 1.0D0))
570 : RNUM = RNUM + ((ALP * GAMMA * PINF) / (GAMMA - 1.0D0))
571 : C
572 : PINF = RNUM / DEN
573 : GAMMA = 1.0D0 / DEN
574 : GAMMA = GAMMA + 1.0D0
575 : PINF = PINF / GAMMA
576 : C
577 : PINFD = PINF
578 : GAMD = GAMMA
579 : C
580 : C A priori, there is no vacuum...
581 : C
582 : LOGVAC = .FALSE.
583 : C
584 : C IF (NGCF .EQ. 685) THEN
585 : C WRITE(*,*)
586 : C & PINFG, GAMG,
587 : C & PINFD, GAMD
588 : C ENDIF
589 : C
590 : IF (INDMET .EQ. 1) THEN
591 : C
592 : C Exact Riemann solver
593 : C Computation of the states on the left and on the
594 : C right of the contact discontinuity
595 : C
595 : CALL FSTIFF(IDIM, NPAR,
596 : & PINFG, GAMG,
597 : & PINFD, GAMD,
598 : & PMIN,
599 : & ROG1, PG1, UNG1, UTG1, 0.0D0, WPARG.PROG,
600 : & ROD1, PD1, UND1, UTD1, 0.0D0, WPARD.PROG,
601 : & RHO_B, P_B, U_B, D_B,
602 : & RHO_A, P_A, U_A, D_A,
603 : & 0.0D0,
604 : & UINT,
605 : & IFLU.PROG, CELLT,
606 : & LOGAN, LOGNC, LOGVAC)
607 : IF (LOGAN) THEN
608 : WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
609 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
610 : C
611 : WRITE(*,*)
612 : WRITE(*,*) 'NGCF =', NGCF
613 : WRITE(*,*) 'NLCF =', NLCF
614 : WRITE(*,*) 'Left '
615 : WRITE(*,*) 'Phic, phi'
616 : WRITE(*,*) PHIG1C, PHIG1F
617 : WRITE(*,*) 'Pinf, gamma'
618 : WRITE(*,*) PINFG, GAMG
619 : WRITE(*,*) 'rho, un, ut, pg'
620 : WRITE(*,*) ROG1, UNG1, UTG1, PG1
621 : WRITE(*,*) 'Right'
622 : WRITE(*,*) 'Phic, phi'
623 : WRITE(*,*) PHID1C, PHID1F
624 : WRITE(*,*) 'Pinf, gamma'
625 : WRITE(*,*) PINFD, GAMD
626 : WRITE(*,*) 'rho, un, ut, pg'
627 : WRITE(*,*) ROD1, UND1, UTD1, PD1
628 : C
628 : CALL ERREUR(5)
629 : GOTO 9999
630 : ENDIF
631 : ENDIF
632 : CELLTM = MAX(CELLT, CELLTM)
633 : C
634 : IF ((.NOT. LOGVAC) .AND. (PHIPRO .GT. 0.0D0)) THEN
635 : C
636 : C************* Ecriture du residu
637 : C
638 : C FLUX(2) = RO Un RO Un
639 : C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
640 : C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
641 : C FLUX(5) = RO Un Et RO Un Et
642 : C Mass fractions... (5 + IESP)...
643 : C Volume fractions in non conservative form...
644 : C
645 : FLUCEL = IFLU.PROG(1) * SURF
646 : C write(*,*) 'flux', iflu.prog(1)
647 : C write(*,*) 'nlceg, nlced', nlceg, nlced
648 : C write(*,*) 'surf', surf
649 : MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
650 : & (FLUCEL / VOLG)
651 : IF (NLCED .NE. NLCEG)
652 : & MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
653 : & (FLUCEL / VOLD)
654 : FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
655 : MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
656 : & (FLUCEL / VOLG)
657 : IF (NLCED .NE. NLCEG)
658 : & MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
659 : & (FLUCEL / VOLD)
660 : FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
661 : MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
662 : & (FLUCEL / VOLG)
663 : IF (NLCED .NE. NLCEG)
664 : & MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
665 : & (FLUCEL / VOLD)
666 : FLUCEL = IFLU.PROG(4) * SURF
667 : MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
668 : & (FLUCEL / VOLG)
669 : IF (NLCED .NE. NLCEG)
670 : & MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
671 : & (FLUCEL / VOLD)
672 : C
673 : C************* rho Y in conservative form
674 : C Components 5 + 1 -> 5 + NESP
675 : C
676 : DO IESP = 1, NESP, 1
677 :
678 : FLUCEL = IFLU.PROG(5 + IESP) * SURF
679 : MPORES.VPOCHA(NLCEG,5 + IESP) =
680 : & MPORES.VPOCHA(NLCEG,5 + IESP) -
681 : & (FLUCEL / VOLG)
682 : IF (NLCED .NE. NLCEG)
683 : & MPORES.VPOCHA(NLCED,5 + IESP) =
684 : & MPORES.VPOCHA(NLCED,5 + IESP) +
685 : & (FLUCEL / VOLD)
686 : ENDDO
687 : C************* Alpha in non conservative form
688 : C Components 5 + NESP -> 5 + (2 * NESP)
689 : IF (NESP .GE. 1) THEN
690 : IF (UINT .GT. 0.0D0) THEN
691 : FRAMAS = ALPHAG
692 : ELSE
693 : FRAMAS = ALPHAD
694 : ENDIF
695 : DO IESP = 1, NESP, 1
696 : C UINT > 0
697 : C FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
698 : C UINT < 0
699 : C FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
700 : FLUCEL = UINT * (FRAMAS.YET(IESP) -
701 : & ALPHCG.YET(IESP)) * SURF
702 : C UINT > 0
703 : C FLUCED = UINT * (PHIG1F - PHID1C) * SURF
704 : C UINT < 0
705 : C FLUCED = UINT * (PHID1F - PHID1C) * SURF
706 : FLUCED = UINT * (FRAMAS.YET(IESP) -
707 : & ALPHCD.YET(IESP)) * SURF
708 : C
709 : MPORES.VPOCHA(NLCEG,5 + NESP + IESP) =
710 : & MPORES.VPOCHA(NLCEG,5 + NESP + IESP) -
711 : & (FLUCEL / VOLG)
712 : IF (NLCED .NE. NLCEG)
713 : & MPORES.VPOCHA(NLCED,5 + NESP + IESP) =
714 : & MPORES.VPOCHA(NLCED,5 + NESP + IESP) +
715 : & (FLUCED / VOLD)
716 : ENDDO
717 : ENDIF
718 : ENDIF
719 : C
720 : C Level set
721 : C
722 : IF (LOGCON) THEN
723 : C
724 : C************* Conservative, variable rho phi...
725 : C
726 : FLUCEL = IFLU.PROG(5) * SURF
727 : C write(*,*) 'flux', iflu.prog(1)
728 : C write(*,*) 'nlceg, nlced', nlceg, nlced
729 : C write(*,*) 'surf', surf
730 : MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
731 : & (FLUCEL / VOLG)
732 : IF (NLCED .NE. NLCEG)
733 : & MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
734 : & (FLUCEL / VOLD)
735 : ELSE
736 : C
737 : C************ Non conservative, variable phi
738 : C
739 : IF (UINT .GT. 0.0D0) THEN
740 : FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
741 : FLUCED = UINT * (PHIG1F - PHID1C) * SURF
742 : ELSE
743 : FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
744 : FLUCED = UINT * (PHID1F - PHID1C) * SURF
745 : ENDIF
746 : C
747 : MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
748 : & (FLUCEL / VOLG)
749 : IF (NLCED .NE. NLCEG)
750 : & MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
751 : & (FLUCED / VOLD)
752 : ENDIF
753 : C
754 : IF ((PHIPRO .LE. 0.0D0) .OR. LOGVAC) THEN
755 : * write(*,*) 'Son qui !!!'
756 : * stop
757 : C
758 : C Left flux for everything but not for the level-set
759 : C
760 : ROD1G = RHO_B
761 : PD1G = P_B
762 : UND1G = U_B
763 : IF (INDMET .EQ. 1) THEN
764 : C
765 : C Exact Riemann solver
766 : C
766 : CALL FSTIFF(IDIM, NPAR,
767 : & PINFG, GAMG,
768 : & PINFG, GAMG,
769 : & PMIN,
770 : & ROG1, PG1, UNG1, UTG1, 0.0D0, WPARG.PROG,
771 : & ROD1G, PD1G, UND1G, UTG1, 0.0D0, WPARG.PROG,
772 : & RHO_B1, P_B1, U_B1, D_B1,
773 : & RHO_A1, P_A1, U_A1, D_A1,
774 : & 0.0D0,
775 : & UINT,
776 : & IFLU.PROG, CELLT,
777 : & LOGAN, LOGNC, LOGVAC)
778 : C IF (LOGVAC) THEN
779 : C WRITE(IOIMP,*) 'Vacuum of vacuum ???'
780 : C WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
781 : C WRITE(IOIMP,*) 'ANOMALY DETECTED'
781 : C CALL ERREUR(251)
782 : C GOTO 9999
783 : C ENDIF
784 : IF (LOGAN) THEN
785 : WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
786 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
787 : C
788 : WRITE(*,*)
789 : WRITE(*,*) 'NGCF =', NGCF
790 : WRITE(*,*) 'Left '
791 : WRITE(*,*) 'Phic, phi'
792 : WRITE(*,*) PHIG1C, PHIG1F
793 : WRITE(*,*) 'rho, un, ut, pg'
794 : WRITE(*,*) ROG1, UNG1, UTG1, PG1
795 : WRITE(*,*) 'Right'
796 : WRITE(*,*) 'Phic, phi'
797 : WRITE(*,*) PHID1C, PHID1F
798 : WRITE(*,*) 'rho, un, ut, pg'
799 : WRITE(*,*) ROD1, UND1, UTD1, PD1
800 : C
800 : CALL ERREUR(5)
801 : GOTO 9999
802 : ENDIF
803 : ENDIF
804 : CELLTM = MAX(CELLT, CELLTM)
805 : C
806 : C********** Ecriture du residu... gauche only...
807 : C
808 : C FLUX(2) = RO Un RO Un
809 : C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
810 : C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
811 : C FLUX(5) = RO Un Et RO Un Et
812 : C
813 : FLUCEL = IFLU.PROG(1) * SURF
814 : C write(*,*) 'flux', iflu.prog(1)
815 : C write(*,*) 'nlceg, nlced', nlceg, nlced
816 : C write(*,*) 'surf', surf
817 : MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
818 : & (FLUCEL / VOLG)
819 : FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
820 : MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
821 : & (FLUCEL / VOLG)
822 : FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
823 : MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
824 : & (FLUCEL / VOLG)
825 : FLUCEL = IFLU.PROG(4) * SURF
826 : MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
827 : & (FLUCEL / VOLG)
828 : C
829 : C Level set already done...
830 : C Level set treated as the other variables (like in the
831 : C following comment... does not work)
832 : C
833 : C FLUCEL = IFLU.PROG(5) * SURF
834 : C MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
835 : C & (FLUCEL / VOLG)
836 : C
837 : C
838 : C************* rho Y in conservative form
839 : C Components 5 + 1 -> 5 + NESP
840 : C
841 : DO IESP = 1, NESP, 1
842 : FLUCEL = IFLU.PROG(5 + IESP) * SURF
843 : MPORES.VPOCHA(NLCEG,5 + IESP) =
844 : & MPORES.VPOCHA(NLCEG,5 + IESP) -
845 : & (FLUCEL / VOLG)
846 : ENDDO
847 : C************* Alpha in non conservative form
848 : C Components 5 + NESP -> 5 + (2 * NESP)
849 : IF (NESP .GE. 1) THEN
850 : IF (UINT .GT. 0.0D0) THEN
851 : FRAMAS = ALPHAG
852 : ELSE
853 : FRAMAS = ALPHAD
854 : ENDIF
855 : DO IESP = 1, NESP, 1
856 : C UINT > 0
857 : C FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
858 : C UINT < 0
859 : C FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
860 : FLUCEL = UINT * (FRAMAS.YET(IESP) -
861 : & ALPHCG.YET(IESP)) * SURF
862 : C
863 : MPORES.VPOCHA(NLCEG,5 + NESP + IESP) =
864 : & MPORES.VPOCHA(NLCEG,5 + NESP + IESP) -
865 : & (FLUCEL / VOLG)
866 : ENDDO
867 : ENDIF
868 : C
869 : C********** Calcul du pas du temps (CFL)
870 : C
871 : C********* a) etat a l'interface
872 : C
873 : DIAMG = MPOVDI.VPOCHA(NLCEG,1)
874 : DIAMD = MPOVDI.VPOCHA(NLCED,1)
875 : DIAM = MIN(DIAMG,DIAMD)
876 : CELL = CELLTM/DIAM
877 : IF(CELL .GT. UNSDT)THEN
878 : UNSDT = CELL
879 : ENDIF
880 : C
881 : C Right flux for everything but not for the level-set
882 : C
883 : ROG1G = RHO_A
884 : PG1G = P_A
885 : UNG1G = U_A
886 : C
887 : IF (NLCED .NE. NLCEG) THEN
888 : IF (INDMET .EQ. 1) THEN
889 : C
890 : C Exact Riemann solver
891 : C
891 : CALL FSTIFF(IDIM, NPAR,
892 : & PINFD, GAMD,
893 : & PINFD, GAMD,
894 : & PMIN,
895 : & ROG1G, PG1G, UNG1G, UTD1, 0.0D0, WPARD.PROG,
896 : & ROD1, PD1, UND1, UTD1, 0.0D0, WPARD.PROG,
897 : & RHO_B1, P_B1, U_B1, D_B1,
898 : & RHO_A1, P_A1, U_A1, D_A1,
899 : & 0.0D0,
900 : & UINT,
901 : & IFLU.PROG, CELLT,
902 : & LOGAN, LOGNC, LOGVAC)
903 : C
904 : C IF (LOGVAC) THEN
905 : C WRITE(IOIMP,*) 'Vacuum of vacuum ???'
906 : C WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
907 : C WRITE(IOIMP,*) 'ANOMALY DETECTED'
907 : C CALL ERREUR(251)
908 : C GOTO 9999
909 : C ENDIF
910 : C
911 : IF (LOGAN) THEN
912 : WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
913 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
913 : CALL ERREUR(5)
914 : GOTO 9999
915 : ENDIF
916 : C
917 : ENDIF
918 : CELLTM = MAX(CELLT, CELLTM)
919 :
920 : C
921 : C************* Ecriture du flux a droite...
922 : C
923 : C FLUX(2) = RO Un RO Un
924 : C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
925 : C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
926 : C FLUX(5) = RO Un Et RO Un Et
927 : C
928 : FLUCEL = IFLU.PROG(1) * SURF
929 : C write(*,*) 'flux', iflu.prog(1)
930 : C write(*,*) 'nlceg, nlced', nlceg, nlced
931 : C write(*,*) 'surf', surf
932 : MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
933 : & (FLUCEL / VOLD)
934 : FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
935 :
936 : MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
937 : & (FLUCEL / VOLD)
938 : FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
939 : MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
940 : & (FLUCEL / VOLD)
941 : FLUCEL = IFLU.PROG(4) * SURF
942 : MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
943 : & (FLUCEL / VOLD)
944 : C
945 : C Level set already done...
946 : C Level set treated as the other variables (like in the
947 : C following comment... does not work)
948 : C
949 : C FLUCEL = IFLU.PROG(5) * SURF
950 : C MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
951 : C & (FLUCEL / VOLD)
952 : C
953 : C
954 : C**************** rho Y in conservative form
955 : C Components 5 + 1 -> 5 + NESP
956 : C
957 : DO IESP = 1, NESP, 1
958 :
959 : FLUCEL = IFLU.PROG(5 + IESP) * SURF
960 : MPORES.VPOCHA(NLCED,5 + IESP) =
961 : & MPORES.VPOCHA(NLCED,5 + IESP) +
962 : & (FLUCEL / VOLD)
963 : ENDDO
964 : C************* Alpha in non conservative form
965 : C Components 5 + NESP -> 5 + (2 * NESP)
966 : IF (NESP .GE. 1) THEN
967 : IF (UINT .GT. 0.0D0) THEN
968 : FRAMAS = ALPHAG
969 : ELSE
970 : FRAMAS = ALPHAD
971 : ENDIF
972 : DO IESP = 1, NESP, 1
973 : FLUCED = UINT * (FRAMAS.YET(IESP) -
974 : & ALPHCD.YET(IESP)) * SURF
975 : C
976 : MPORES.VPOCHA(NLCED,5 + NESP + IESP) =
977 : & MPORES.VPOCHA(NLCED,5 + NESP + IESP) +
978 : & (FLUCED / VOLD)
979 : ENDDO
980 : ENDIF
981 : ENDIF
982 : ENDIF
983 : DIAMG = MPOVDI.VPOCHA(NLCEG,1)
984 : DIAMD = MPOVDI.VPOCHA(NLCED,1)
985 : DIAM = MIN(DIAMG,DIAMD)
986 : CELL = CELLTM/DIAM
987 : IF(CELL .GT. UNSDT)THEN
988 : UNSDT = CELL
989 : ENDIF
990 : ENDIF
991 : C
992 : C**** Fin boucle sur FACEL
993 : C
994 : ENDDO
995 : C
996 : C**** Pas du temps (condition de non interaction en 1D)
997 : C
998 : DT = 0.5D0 / UNSDT
999 : C
1000 : C**** Desactivation des segments
1001 : C
1002 : SEGSUP MLELIM
1003 : SEGSUP MLENT1
1004 : SEGSUP MLENT2
1005 : C
1006 : SEGDES MPOVSU
1007 : SEGDES MPOVDI
1008 : SEGDES MPOVVO
1009 : SEGDES MPORES
1010 : SEGDES IPT2
1011 : C
1012 : C Warning. MLRMFR = PGAS . 'MASSFRA'
1013 : C MLRCHE = PGAS . 'CHEMCOEF'
1014 : C
1015 : SEGSUP WPARG
1016 : SEGSUP WPARD
1017 : C
1018 : SEGSUP IFLU
1019 : C
1020 : SEGDES MPPHI1
1021 : SEGDES MELPH1
1022 : SEGDES MELRO1
1023 : SEGDES MELP1
1024 : SEGDES MELUN1
1025 : SEGDES MELUT1
1026 : SEGDES MELVNX
1027 : SEGDES MELVNY
1028 : SEGDES MELT1X
1029 : SEGDES MELT1Y
1030 : C
1031 : SEGDES MLRMGA
1032 : SEGDES MLRPGA
1033 : SEGDES MLRMPI
1034 : SEGDES MLRPPI
1035 : C
1036 : IF (NESP .GE. 1)THEN
1037 : C
1038 : SEGDES MPALPH
1039 : C
1040 : SEGSUP FRAMAG
1041 : SEGSUP FRAMAD
1042 : SEGSUP ALPHAG
1043 : SEGSUP ALPHAD
1044 : SEGSUP ALPHCG
1045 : SEGSUP ALPHCD
1046 : C
1047 : DO IESP = 1, NESP
1048 : MELVA1 = MCHAMY.IELVAL(IESP)
1049 : SEGDES MELVA1
1050 : ENDDO
1051 : SEGDES MCHAMY
1052 : C
1053 : DO IESP = 1, NESP
1054 : MELVA1 = MCHAMA.IELVAL(IESP)
1055 : SEGDES MELVA1
1056 : ENDDO
1057 : SEGDES MCHAMA
1058 : ENDIF
1059 : C
1060 : 9999 CONTINUE
1061 : RETURN
1062 : END
1063 : CC
1064 :
1065 :
1066 :
© Cast3M 2003 - All rights reserved.
Disclaimer