1 : C KODFL1 SOURCE BECC 11/06/08 21:15:36 7000
2 : SUBROUTINE KODFL1(INDK0,NESP,NORD,TMAX,RUNIV,PROPHY,
3 : & MLRCHE,MLRMFR,
4 : & INDMET,
5 : & ICHPA1, ICHPA2,
6 : & IALF1,IROF1,IVITF1,IPF1,
7 : & IALF2,IROF2,IVITF2,IPF2,
8 : & IUINF1, IUINF2,
9 : & K0, ICHPK0,IGRALP,EPS,
10 : & ICHPSU,ICHPDI,ICHPVO,
11 : & MELEMC,MELEMF,MELEFE,MELLIM,
12 : & ICHRES,
13 : & DT,SURFL,
14 : & LOGNC,LOGAN,MESERR)
15 : C************************************************************************
16 : C
17 : C PROJET : CASTEM 2000
18 : C
19 : C NOM : KODFL1
20 : C
21 : C DESCRIPTION : Voir KONV13
22 : C
23 : C Cas deux dimensions, DEM
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 0) INDK0 = 1 => K0 = constant
34 : C INDK0 = 2 => K0 = variable
35 : C
36 : C 1) Parameters concerning some gas properties
37 : C
38 : C NESP : nombre d'especes dans le melange.
39 : C
40 : C NORD : ordre des polynoms du cv_i
41 : C
42 : C TMAX : maximum temperature for cv expansion
43 : C
44 : C RUNIV : universal constant for gases
45 : C
46 : C PROPHY : thermodynamic properties for the gases
47 : C
48 : C MLRCHE : LISTREEL with the coefS involved in the chemical
49 : C reaction
50 : C
51 : C MLRMFR : LISTREEL with the mass fractions before or after
52 : C the chemical reaction
53 : C
54 : C 2) Numerical scheme
55 : C
56 : C INDMET : 1 VLH (non-reactive) + SS (reactive)
57 : C
58 : C INDMET : 2 SS (non-reactive) + SS (reactive)
59 : C
60 : C INDMET : 3 AUSM+up (non-reactive) + SS (reactive)
61 : C
62 : C 3) Pointers for the CHPOINTs/MCHAML des variables
63 : C
64 : C ICHPA1 : chpoint CENTRE contenant alpha1
65 : C ICHPA2 : chpoint CENTRE contenant alpha2
66 : C
67 : C IALF1 : MCHAML sur "FACEL" contenant alpha1
68 : C ("gauche" et "droite");
69 : C
70 : C IROF1 : MCHAML sur "FACEL" contenant la masse volumique 1
71 : C ("gauche" et "droite");
72 : C
73 : C IVITF1 : MCHAML sur "FACEL" contenant la vitesse 1 dans le repaire
74 : C local (n,t) et les cosinus directeurs des repaire local;
75 : C
76 : C IPF1 : MCHAML sur "FACEL" contenant la pression 1;
77 : C
78 : C IALF2 : MCHAML sur "FACEL" contenant alpha1
79 : C ("gauche" et "droite");
80 : C
81 : C IROF2 : MCHAML sur "FACEL" contenant la masse volumique 1
82 : C ("gauche" et "droite");
83 : C
84 : C IVITF2 : MCHAML sur "FACEL" contenant la vitesse 1 dans le repaire
85 : C local (n,t) et les cosinus directeurs des repaire local;
86 : C
87 : C IPF2 : MCHAML sur "FACEL" contenant la pression 1;
88 : C
89 : C IUINF1 : CHPOINT, one component, "cut-off velocity" 1
90 : C 0 if no Bas MACH
91 : C
92 : C IUINF2 : CHPOINT, one component, "cut-off velocity" 2
93 : C
94 : C
95 : C 4) Others
96 : C
97 : C K0 or
98 : C ICHPK0 : fundamental flame speed
99 : C
100 : C IGRALP : CHPOINT "FACE" contenant grad(alpha)/alpha
101 : C
102 : C EPS : "zero"
103 : C
104 : C ICHPSU : CHPOINT "FACE" contenant la surface des faces
105 : C
106 : C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
107 : C de chaque element
108 : C
109 : C ICHPVO : CHPOINT "CENTRE" contenant le volume
110 : C de chaque element
111 : C
112 : C 5) Pointeurs de MELEME de la table DOMAINE
113 : C
114 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
115 : C
116 : C MELEMF : MELEME 'FACE' du SPG des FACES
117 : C
118 : C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
119 : C
120 : C MELLIM : MAILLAGE where the flux (or residu) will not be found
121 : C
122 : C SORTIES (il faudrait dire E/S)
123 : C
124 : C ICHRES : pointeurs de CHPOINTs "FACE" des flux aux interfaces:
125 : C
126 : C DT : pas de temps pour le respect de la CFL-like condition
127 : C DT < DIAMEL /2 /max(Lambda_i)
128 : C En maillage regulier cette condition garantie la
129 : C non-interaction des ondes
130 : C
131 : C LOGNC : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
132 : C dans pour la solution du probleme Riemann n'a pas bien
133 : C marchéee
134 : C
135 : C LOGAN : (LOGICAL): si .TRUE. une anomalie à été detectée
136 : C
137 : C MESERR : pour l'ecriture des messages d'erreurs
138 : C
139 : C************************************************************************
140 : C
141 : C HISTORIQUE (Anomalies et modifications éventuelles)
142 : C
143 : C 07.12.09 Creation
144 : C 21.12.2010 Estension au 3D
145 : C
146 : C************************************************************************
147 : C
148 : C
149 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0
150 : C Y \in (0,1)
151 : C Si non il faut le faire!!!
152 : C
153 : C************************************************************************
154 : C
155 : C
156 : IMPLICIT INTEGER(I-N)
157 : INTEGER I1, I2, ICOMP
158 : & , NESP,NORD ,INDMET
159 : & , ICHPA1, ICHPA2
160 : & , IALF1, IROF1, IVITF1, IPF1
161 : & , IALF2, IROF2, IVITF2, IPF2
162 : & , IUINF1, IUINF2
163 : & , IGRALP, ICHPK0
164 : & , ICHPSU, ICHPDI, ICHPVO, MELEMC, MELEMF, MELEFE, MELLIM
165 : & , IGEOMC, IGEOMF
166 : & , ICHRES
167 : & , NFAC
168 : & , NLCF, NGCEG, NGCED, NLCEG, NLCED
169 : & , NGCF, NLCF1, SPG1, SPG2, NLFL
170 : & , NLHS
171 : REAL*8 TMAX, RUNIV, K0, K0G, K0D, EPS
172 : REAL*8 DT, UNSDT, CELLT, CELLTM
173 : & , AG1C, ALPG1, ROG1, UNG1, UTG1, UVG1, PG1
174 : & , AG2C, ALPG2, ROG2, UNG2, UTG2, UVG2, PG2
175 : & , AD1C, ALPD1, ROD1, UND1, UTD1, UVD1, PD1
176 : & , AD2C, ALPD2, ROD2, UND2, UTD2, UVD2, PD2
177 : & , SURF,CNX, CNY, CNZ, CTX , CTY, CTZ, CVX, CVY, CVZ
178 : & , V_INF
179 : & , GRALX, GRALY, GRALZ
180 : & , CNNA, CTNA, CVNA
181 : & , CELL, DIAMG, DIAMD, DIAM, VOLG, VOLD
182 : & , FLUCEL, SCON1, SG1, SD1, SCON2, SG2, SD2
183 : & , S1_11, S2_22, S1_12, S1_21, S2_12, S2_21
184 : & , S1_21l, S1_12l, S2_21l, S2_12l
185 : & , S1_21r, S1_12r, S2_21r, S2_12r
186 : & , S1_21e, S1_12e, S2_21e, S2_12e
187 : & , RESG(12), RESD(12)
188 : & , SURFL, SURF0
189 : LOGICAL LOGNC, LOGAN
190 : CHARACTER*(40) MESERR
191 : CHARACTER*(8) TYPE
192 : C
193 : C**** Segment des proprietes du gaz
194 : C
195 : SEGMENT PROPHY
196 : REAL*8 XTAB(N1,N2)
197 : C REAL*8 ACV(NORDP1,NESP+1), W(NESP+1), H0K(NESP+1)
198 : C & ,ACVTOG(NORDP1), ACVTOD(NORDP1)
199 : ENDSEGMENT
200 : C N1 = NORD + 3
201 : C N2 = NESP
202 : C
203 : C**** LES INCLUDES
204 : C
205 : -INC CCOPTIO
206 : -INC SMCOORD
207 : -INC SMELEME
208 : -INC SMCHAML
209 : C Normal and tangent
210 : POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL,
211 : & MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL,
212 : & MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL
213 : POINTEUR MELUN1.MELVAL, MELUT1.MELVAL, MELUV1.MELVAL
214 : POINTEUR MELAL1.MELVAL, MELRO1.MELVAL, MELP1.MELVAL
215 : POINTEUR MELUN2.MELVAL, MELUT2.MELVAL, MELUV2.MELVAL
216 : POINTEUR MELAL2.MELVAL, MELRO2.MELVAL, MELP2.MELVAL
217 : -INC SMCHPOI
218 : POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL, MPOVVO.MPOVAL
219 : & , MPORES.MPOVAL, MPOALP.MPOVAL
220 : & , MPALP1.MPOVAL, MPALP2.MPOVAL
221 : & , MPUIN2.MPOVAL, MPUIN1.MPOVAL
222 : & , MPOK0.MPOVAL
223 : -INC SMLMOTS
224 : -INC SMLENTI
225 : POINTEUR MLELIM.MLENTI
226 : INTEGER JG
227 : -INC SMLREEL
228 : POINTEUR MLRMFR.MLREEL, MLRCHE.MLREEL
229 : C
230 : C**** Segment for the subroutine prithe
231 : C
232 : POINTEUR IWST_L.MLREEL, IWST_R.MLREEL, IWWORK.MLREEL
233 : C 1:(4+IDIM+NESP)
234 : POINTEUR IY.MLREEL
235 : C 1:NESP
236 : POINTEUR ICVL.MLREEL, ICVR.MLREEL, ICVRSS.MLREEL, ICVWOR.MLREEL
237 : C 0:NORD
238 : C
239 : POINTEUR IFLU.MLREEL, IFLULA.MLREEL
240 : C
241 : C**** Surface de flamme totale
242 : C
243 : SURFL = 0.0D0
244 : C
245 : C**** KRIPAD pour la correspondance global/local des conditions limits
246 : C
246 : CALL KRIPAD(MELLIM,MLELIM)
247 : C SEGINI MLELIM
248 : C
249 : C**** Initialisation des MCHAMLs
250 : C
251 : C Alpha
252 : C
253 : MCHEL1 = IALF1
254 : SEGACT MCHEL1
255 : MCHAM1 = MCHEL1.ICHAML(1)
256 : SEGACT MCHAM1
257 : MELAL1 = MCHAM1.IELVAL(1)
258 : SEGDES MCHEL1
259 : SEGDES MCHAM1
260 : C
261 : C Masse volumique
262 : C
263 : MCHEL1 = IROF1
264 : SEGACT MCHEL1
265 : MCHAM1 = MCHEL1.ICHAML(1)
266 : SEGACT MCHAM1
267 : MELRO1 = MCHAM1.IELVAL(1)
268 : SEGDES MCHEL1
269 : SEGDES MCHAM1
270 : C
271 : C Pression
272 : C
273 : MCHEL1 = IPF1
274 : SEGACT MCHEL1
275 : MCHAM1 = MCHEL1.ICHAML(1)
276 : SEGACT MCHAM1
277 : MELP1 = MCHAM1.IELVAL(1)
278 : SEGDES MCHEL1
279 : SEGDES MCHAM1
280 : C
281 : C Vitesse et cosinus directeurs du repere (n,t)
282 : C
283 : MCHEL1 = IVITF1
284 : SEGACT MCHEL1
285 : C
286 : C La vitesse a comme SPG MELEFE
287 : C Le cosinus directeurs ont comme SPG MELEMF
288 : C
289 : C MCHAM1 -> Cosinus directeurs
290 : C MCHAM2 -> Vitesse
291 : C
292 : SPG1 = MCHEL1.IMACHE(1)
293 : SPG2 = MCHEL1.IMACHE(2)
294 : IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
295 : MCHAM1 = MCHEL1.ICHAML(1)
296 : MCHAM2 = MCHEL1.ICHAML(2)
297 : ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
298 : MCHAM1 = MCHEL1.ICHAML(2)
299 : MCHAM2 = MCHEL1.ICHAML(1)
300 : ELSE
301 : LOGAN = .TRUE.
302 : GOTO 9999
303 : ENDIF
304 : SEGACT MCHAM1
305 : SEGACT MCHAM2
306 : IF (IDIM .EQ. 2) THEN
307 : MELVNX = MCHAM1.IELVAL(1)
308 : MELVNY = MCHAM1.IELVAL(2)
309 : MELT1X = MCHAM1.IELVAL(3)
310 : MELT1Y = MCHAM1.IELVAL(4)
311 : MELUN1 = MCHAM2.IELVAL(1)
312 : MELUT1 = MCHAM2.IELVAL(2)
313 : ELSE
314 : MELVNX = MCHAM1.IELVAL(1)
315 : MELVNY = MCHAM1.IELVAL(2)
316 : MELVNZ = MCHAM1.IELVAL(3)
317 : MELT1X = MCHAM1.IELVAL(4)
318 : MELT1Y = MCHAM1.IELVAL(5)
319 : MELT1Z = MCHAM1.IELVAL(6)
320 : MELT2X = MCHAM1.IELVAL(7)
321 : MELT2Y = MCHAM1.IELVAL(8)
322 : MELT2Z = MCHAM1.IELVAL(9)
323 : MELUN1 = MCHAM2.IELVAL(1)
324 : MELUT1 = MCHAM2.IELVAL(2)
325 : MELUV1 = MCHAM2.IELVAL(3)
326 : ENDIF
327 : SEGDES MCHAM1
328 : SEGDES MCHAM2
329 : SEGDES MCHEL1
330 : C
331 : C Alpha
332 : C
333 : MCHEL1 = IALF2
334 : SEGACT MCHEL1
335 : MCHAM1 = MCHEL1.ICHAML(1)
336 : SEGACT MCHAM1
337 : MELAL2 = MCHAM1.IELVAL(1)
338 : SEGDES MCHEL1
339 : SEGDES MCHAM1
340 : C
341 : C Masse volumique
342 : C
343 : MCHEL1 = IROF2
344 : SEGACT MCHEL1
345 : MCHAM1 = MCHEL1.ICHAML(1)
346 : SEGACT MCHAM1
347 : MELRO2 = MCHAM1.IELVAL(1)
348 : SEGDES MCHEL1
349 : SEGDES MCHAM1
350 : C
351 : C Pression
352 : C
353 : MCHEL1 = IPF2
354 : SEGACT MCHEL1
355 : MCHAM1 = MCHEL1.ICHAML(1)
356 : SEGACT MCHAM1
357 : MELP2 = MCHAM1.IELVAL(1)
358 : SEGDES MCHEL1
359 : SEGDES MCHAM1
360 : C
361 : C Vitesse et cosinus directeurs du repere (n,t)
362 : C
363 : MCHEL1 = IVITF2
364 : SEGACT MCHEL1
365 : C
366 : C La vitesse a comme SPG MELEFE
367 : C Le cosinus directeurs ont comme SPG MELEMF
368 : C
369 : C MCHAM1 -> Cosinus directeurs
370 : C MCHAM2 -> Vitesse
371 : C
372 : SPG1 = MCHEL1.IMACHE(1)
373 : SPG2 = MCHEL1.IMACHE(2)
374 : IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
375 : MCHAM1 = MCHEL1.ICHAML(1)
376 : MCHAM2 = MCHEL1.ICHAML(2)
377 : ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
378 : MCHAM1 = MCHEL1.ICHAML(2)
379 : MCHAM2 = MCHEL1.ICHAML(1)
380 : ELSE
381 : LOGAN = .TRUE.
382 : GOTO 9999
383 : ENDIF
384 : C SEGACT MCHAM1
385 : C MELVNX = MCHAM1.IELVAL(1)
386 : C MELVNY = MCHAM1.IELVAL(2)
387 : C MELT1X = MCHAM1.IELVAL(3)
388 : C MELT1Y = MCHAM1.IELVAL(4)
389 : C SEGDES MCHAM1
390 : SEGACT MCHAM2
391 : MELUN2 = MCHAM2.IELVAL(1)
392 : MELUT2 = MCHAM2.IELVAL(2)
393 : IF (IDIM .EQ. 3) MELUV2 = MCHAM2.IELVAL(3)
394 : SEGDES MCHAM2
395 : SEGDES MCHEL1
396 : C write(*,*) 'Son in kodfl1'
397 : C write(*,*) 'Ho letto gli MCHAM'
398 : C
399 : C**** Segment for fcolre
400 : C
401 : C Chemical coefficients in MLRCHE.MLREEL.
402 : SEGACT MLRCHE
403 : NLHS = 0
404 : DO I1 = 1, NESP, 1
405 : IF (MLRCHE.PROG(I1) .GT. 0.0D0) NLHS = NLHS + 1
406 : C write(*,*) MLRCHE.PROG(I1)
407 : ENDDO
408 : c write(*,*) NLHS
409 : C
410 : C Mass fractions
411 : C MLRMFR.PROG(1) = YINI_first
412 : C MLRMFR.PROG(2) = YFIN_first
413 : C MLRMFR.PROG(3) = YFIN_second
414 : C ...
415 : C MLRMFR.PROG(nLHS+1) = YFIN_nLHS
416 : C MLRMFR.PROG(nLHS+2) = YINI_(nLHS+1)
417 : C
418 : SEGACT MLRMFR
419 : C
420 : C cv, molar weicht formation energy in PROPHY.XTAB,
421 : C already active
422 : C
423 : C States
424 : JG = NESP + IDIM + 4
425 : SEGINI IWWORK
426 : SEGINI IWST_R
427 : SEGINI IWST_L
428 : C
429 : C Mass fractions
430 : JG = NESP
431 : SEGINI IY
432 : C
433 : C Coeff for CV
434 : JG = NORD + 1
435 : SEGINI ICVL
436 : SEGINI ICVR
437 : SEGINI ICVRSS
438 : SEGINI ICVWOR
439 : C
440 : C Fluxes
441 : JG = NESP + IDIM + 4
442 : SEGINI IFLU
443 : SEGINI IFLULA
444 : C
445 : C Since the mixture here considered is homogeneous, we first fill
446 : C the part of IWST_L and IWST_R involving the species mass fractions
447 : C and K0
448 : C
449 : IWST_L.PROG(4 + IDIM) = MLRMFR.PROG(2)
450 : IWST_R.PROG(4 + IDIM) = MLRMFR.PROG(2)
451 : DO I2 = 3, NESP, 1
452 : I1 = 2 + IDIM + I2
453 : C DO I1 = 5 + IDIM, 2 + IDIM + NESP
454 : IWST_L.PROG(I1) = MLRMFR.PROG(I2)
455 : IWST_R.PROG(I1) = MLRMFR.PROG(I2)
456 : ENDDO
457 : IWST_L.PROG(3 + IDIM + NESP) = MLRMFR.PROG(1) - MLRMFR.PROG(2)
458 : IWST_R.PROG(3 + IDIM + NESP) = MLRMFR.PROG(1) - MLRMFR.PROG(2)
459 : C
460 : C
461 : C**** Initialisation des MELEMEs
462 : C
463 : C 'CENTRE', 'FACEL'
464 : C
465 : IPT2 = MELEFE
466 : SEGACT IPT2
467 : NFAC = IPT2.NUM(/2)
468 : C
469 : C**** KRIPAD pour la correspondance global/local de centre
470 : C
470 : CALL KRIPAD(MELEMC,MLENT1)
471 : C
472 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
473 : C
474 : C Si i est le numero global d'un noeud de ICEN,
475 : C MLENT1.LECT(i) contient sa position, i.e.
476 : C
477 : C I = numero global du noeud centre
478 : C MLENT1.LECT(i) = numero local du noeud centre
479 : C
480 : C MLENT1 déjà activé, i.e.
481 : C
482 : C SEGINI MLENT1
483 : C
484 : C
485 : C**** KRIPAD pour la correspondance global/local de 'FACE'
486 : C
486 : CALL KRIPAD(MELEMF,MLENT2)
487 : C SEGINI MLENT2
488 : C
489 : C
490 : C**** CHPOINTs de la normale à alpha
491 : C
491 : CALL LICHT(IGRALP,MPOALP,TYPE,IGEOMF)
491 : IF (INDK0 .EQ. 2) CALL LICHT(ICHPK0,MPOK0,TYPE,IGEOMF)
492 : C SEGACT MPOALP
493 : C SEGDES MPOK0
494 : C
495 : C**** CHPOINTs de la table DOMAINE
496 : C
496 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
496 : CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
496 : CALL LICHT(ICHPVO,MPOVVO,TYPE,IGEOMC)
497 : C
498 : C**** LICHT active les MPOVALs en *MOD
499 : C
500 : C i.e.
501 : C
502 : C SEGACT MPOVSU*MOD
503 : C SEGACT MPOVDI*MOD
504 : C SEGACT MPOVVO*MOD
505 : C
506 : C
507 : C**** Les FLUX aux face
508 : C
508 : CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
509 : C SEGACT MPORES*MOD
510 : C
511 : C**** Activation des MCHAMLs
512 : C
513 : SEGACT MELAL1
514 : SEGACT MELRO1
515 : SEGACT MELP1
516 : SEGACT MELUN1
517 : SEGACT MELUT1
518 : SEGACT MELAL2
519 : SEGACT MELRO2
520 : SEGACT MELP2
521 : SEGACT MELUN2
522 : SEGACT MELUT2
523 : SEGACT MELVNX
524 : SEGACT MELVNY
525 : SEGACT MELT1X
526 : SEGACT MELT1Y
527 : IF (IDIM .EQ. 3) THEN
528 : SEGACT MELUV1
529 : SEGACT MELUV2
530 : SEGACT MELVNZ
531 : SEGACT MELT1Z
532 : SEGACT MELT2X
533 : SEGACT MELT2Y
534 : SEGACT MELT2Z
535 : ENDIF
536 : C
537 : C**** Lecture des MPOVALs ICHPA1, ICHPA2
538 : C
538 : CALL LICHT(ICHPA1, MPALP1, TYPE, IGEOMC)
538 : CALL LICHT(ICHPA2, MPALP2, TYPE, IGEOMC)
539 : C SEGACT MPALP1
540 : C SEGACT MPALP2
541 : C
542 : C**** Bas Mach
543 : C
544 : IF(IUINF1 .NE. 0)THEN
544 : CALL LICHT(IUINF2,MPUIN2,TYPE,IGEOMC)
544 : CALL LICHT(IUINF1,MPUIN1,TYPE,IGEOMC)
545 : C SEGACT MPUIN2*MOD
546 : C SEGACT MPUIN1*MOD
547 : ELSE
548 : MPUIN2=0
549 : MPUIN1=0
550 : ENDIF
551 : C
552 : C**** Initialisation de 1/DT et d'autres variables
553 : C
554 : UNSDT = 0.0D0
555 : GRALZ = 0.0D0
556 : CNZ = 0.0D0
557 : CTZ = 0.0D0
558 : CVX = 0.0D0
559 : CVY = 0.0D0
560 : CVZ = 0.0D0
561 : UVG1 = 0.0D0
562 : UVD1 = 0.0D0
563 : UVG2 = 0.0D0
564 : UVD2 = 0.0D0
565 : C
566 : C**** BOUCLE SUR FACEL pour le calcul du FLUX/RESIDU
567 : C
568 : DO NLCF = 1, NFAC
569 : C
570 : C******* NLCF = numero local du centre de facel
571 : C NGCF = numero global du centre de facel
572 : C NLCF1 = numero local du centre de face
573 : C NGCEG = numero global du centre ELT "gauche"
574 : C NLCEG = numero local du centre ELT "gauche"
575 : C NGCED = numero global du centre ELT "droite"
576 : C NLCED = numero local du centre ELT "droite"
577 : C
578 : NGCEG = IPT2.NUM(1,NLCF)
579 : NGCED = IPT2.NUM(3,NLCF)
580 : NGCF = IPT2.NUM(2,NLCF)
581 : NLCF1 = MLENT2.LECT(NGCF)
582 : NLCEG = MLENT1.LECT(NGCEG)
583 : NLCED = MLENT1.LECT(NGCED)
584 : NLFL = MLELIM.LECT(NGCF)
585 : C
586 : VOLG = MPOVVO.VPOCHA(NLCEG,1)
587 : VOLD = MPOVVO.VPOCHA(NLCED,1)
588 : SURF0 = MPOVSU.VPOCHA(NLCF,1)
589 : C
590 : C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
591 : C
592 : IF(NLCF .NE. NLCF1)THEN
593 : MESERR = 'Il ne faut pas jouer avec la console. '
594 : LOGAN = .TRUE.
595 : GOTO 9999
596 : ENDIF
597 : IF(NLFL .EQ. 0)THEN
598 : CELLTM = 0.0D0
599 : C
600 : C********** Recuperation des Etats "gauche" et "droite"
601 : C
602 : AG1C = MPALP1.VPOCHA(NLCEG,1)
603 : ALPG1 = MELAL1.VELCHE(1,NLCF)
604 : ROG1 = MELRO1.VELCHE(1,NLCF)
605 : UNG1 = MELUN1.VELCHE(1,NLCF)
606 : UTG1 = MELUT1.VELCHE(1,NLCF)
607 : PG1 = MELP1.VELCHE(1,NLCF)
608 : AG2C = MPALP2.VPOCHA(NLCEG,1)
609 : ALPG2 = MELAL2.VELCHE(1,NLCF)
610 : ROG2 = MELRO2.VELCHE(1,NLCF)
611 : UNG2 = MELUN2.VELCHE(1,NLCF)
612 : UTG2 = MELUT2.VELCHE(1,NLCF)
613 : PG2 = MELP2.VELCHE(1,NLCF)
614 : C
615 : AD1C = MPALP1.VPOCHA(NLCED,1)
616 : ALPD1 = MELAL1.VELCHE(3,NLCF)
617 : ROD1 = MELRO1.VELCHE(3,NLCF)
618 : UND1 = MELUN1.VELCHE(3,NLCF)
619 : UTD1 = MELUT1.VELCHE(3,NLCF)
620 : PD1 = MELP1.VELCHE(3,NLCF)
621 : AD2C = MPALP2.VPOCHA(NLCED,1)
622 : ALPD2 = MELAL2.VELCHE(3,NLCF)
623 : ROD2 = MELRO2.VELCHE(3,NLCF)
624 : UND2 = MELUN2.VELCHE(3,NLCF)
625 : UTD2 = MELUT2.VELCHE(3,NLCF)
626 : PD2 = MELP2.VELCHE(3,NLCF)
627 : C
628 : CNX = MELVNX.VELCHE(1,NLCF)
629 : CNY = MELVNY.VELCHE(1,NLCF)
630 : CTX = MELT1X.VELCHE(1,NLCF)
631 : CTY = MELT1Y.VELCHE(1,NLCF)
632 : C
633 : GRALX = MPOALP.VPOCHA(NLCF,1)
634 : GRALY = MPOALP.VPOCHA(NLCF,2)
635 : C
636 : C GRALX = CNX
637 : C GRALY = CNY
638 : C
639 : IF (IDIM .EQ. 3) THEN
640 : GRALZ = MPOALP.VPOCHA(NLCF,3)
641 : CNZ = MELVNZ.VELCHE(1,NLCF)
642 : CTZ = MELT1Z.VELCHE(1,NLCF)
643 : CVX = MELT2X.VELCHE(1,NLCF)
644 : CVY = MELT2Y.VELCHE(1,NLCF)
645 : CVZ = MELT2Z.VELCHE(1,NLCF)
646 : UVG1 = MELUV1.VELCHE(1,NLCF)
647 : UVD1 = MELUV1.VELCHE(3,NLCF)
648 : UVG2 = MELUV2.VELCHE(1,NLCF)
649 : UVD2 = MELUV2.VELCHE(3,NLCF)
650 : ENDIF
651 : C
652 : CNNA = (GRALX * CNX) + (GRALY * CNY) + (GRALZ * CNZ)
653 : CTNA = (GRALX * CTX) + (GRALY * CTY) + (GRALZ * CTZ)
654 : CVNA = (GRALX * CVX) + (GRALY * CVY) + (GRALZ * CVZ)
655 : C
656 : C if ((ngcf .eq. 90) .or. (ngcf .eq. 84) .or. (ngcf .eq. 110))
657 : C $ then
658 : C write(*,*) 'ngcf = ', ngcf
659 : C endif
660 : C
661 : C write(*,*) cnna
662 : C
663 : C write(*,*) alpg1, rog1, ung1, utg1, pg1
664 : C write(*,*) alpg2, rog2, ung2, utg2, pg2
665 : C write(*,*) alpd1, rod1, und1, utd1, pd1
666 : C write(*,*) alpd2, rod2, und2, utd2, pd2
667 : C
668 : C
669 : C
670 : C********** Computation of the partition of the surface
671 : C
672 : S1_11 = min (ALPG1, ALPD1)
673 : S2_22 = min (ALPG2, ALPD2)
674 : C S2_22 = 1.0d0 - (max (ALPG1, ALPD1))
675 : S1_21 = max (0.0d0, (AD1C - AG1C))
676 : S1_21r = max (0.0d0, (AD1C - ALPD1))
677 : S1_21e = max (0.0d0, (ALPD1 - ALPG1))
678 : S1_21l = max (0.0d0, (ALPG1 - AG1C))
679 : S2_21 = max (0.0d0, (AG2C - AD2C))
680 : S2_21l = max (0.0d0, (AG2C - ALPG2))
681 : S2_21e = max (0.0d0, (ALPG2 - ALPD2))
682 : S2_21r = max (0.0d0, (ALPD2 - AD2C))
683 : C
684 : C S1_21 = min (S1_21, S2_21)
685 : C S2_21 = S1_21
686 : S1_12 = max (0.0d0, (AG1C - AD1C))
687 : S1_12l = max (0.0d0, (AG1C - ALPG1))
688 : S1_12e = max (0.0d0, (ALPG1 - ALPD1))
689 : S1_12r = max (0.0d0, (ALPD1 - AD1C))
690 : S2_12 = max (0.0d0, (AD2C - AG2C))
691 : S2_12r = max (0.0d0, (AD2C - ALPD2))
692 : S2_12e = max (0.0d0, (ALPD2 - ALPG2))
693 : S2_12l = max (0.0d0, (ALPG2 - AG2C))
694 : C
695 : C write(*,*) 'NGCF =', NGCF
696 : C write(*,*) (S1_11 + S2_22 + S1_12 + S1_21)
697 : C write(*,*) (S1_11 + S2_22 + S2_12 + S2_21)
698 : *
699 : C write(*,*) (S1_11 + S2_22 + S1_12 + S1_21)
700 : C write(*,*) (S1_11 + S2_22 + S2_12 + S2_21)
701 : c CNNA = 1.0D0
702 : c CTNA = 0.0D0
703 : c CVNA = 0.0D0
704 : C
705 : C SURFL = Surface which burns * (CNNA)
706 : *
707 : SURFL = SURFL + ((ABS ((S1_12 + S1_21) * CNNA)) * SURF0)
708 : C
709 : C Flame speed
710 : C
711 : IF (INDK0 .EQ. 1) THEN
712 : K0G = K0
713 : K0D = K0
714 : ELSE
715 : K0G = MPOK0.VPOCHA(NLCEG,1)
716 : K0D = MPOK0.VPOCHA(NLCED,1)
717 : ENDIF
718 : IWST_L.PROG(4 + IDIM + NESP) = K0G
719 : IWST_R.PROG(4 + IDIM + NESP) = K0D
720 : C
721 : C********** Riemann problem 1-1
722 : C
723 : IF (S1_11 .GE. EPS) THEN
724 : IWST_L.PROG(1) = ROG1
725 : IWST_L.PROG(2) = UNG1
726 : C IWST_L.PROG(2) = 0.0D0
727 : IWST_L.PROG(3) = UTG1
728 : C IDIM + 2 = 2 in 2D case
729 : IWST_L.PROG(IDIM + 2)= PG1
730 : IWST_L.PROG(IDIM + 3) = 0.0D0
731 : C
732 : IWST_R.PROG(1) = ROD1
733 : IWST_R.PROG(2) = UND1
734 : C IWST_R.PROG(2) = 0.0D0
735 : IWST_R.PROG(3) = UTD1
736 : IWST_R.PROG(IDIM + 2) = PD1
737 : IWST_R.PROG(IDIM + 3) = 0.0D0
738 : IF (IDIM .EQ. 3) THEN
739 : IWST_L.PROG(IDIM + 1) = UVG1
740 : IWST_R.PROG(IDIM + 1) = UVD1
741 : ENDIF
742 : C
743 : C do i2 = 1 , 4 + NESP + IDIM, 1
744 : C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
745 : C enddo
746 : C
747 : IF (INDMET .EQ. 1) THEN
747 : CALL FVLHRE(IDIM, NESP, NLHS, NORD,
748 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
749 : & TMAX,
750 : & CNNA, CTNA, CVNA,
751 : & IWST_L.PROG,
752 : & IWST_R.PROG,
753 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
754 : & ICVRSS.PROG, ICVWOR.PROG,
755 : & IFLU.PROG, IFLULA.PROG, CELLT,
756 : & LOGAN)
757 : ELSEIF(INDMET .EQ. 2)THEN
757 : CALL FBECRE(IDIM, NESP, NLHS, NORD,
758 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
759 : & TMAX,
760 : & CNNA, CTNA, CVNA,
761 : & IWST_L.PROG,
762 : & IWST_R.PROG,
763 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
764 : & ICVRSS.PROG, ICVWOR.PROG,
765 : & IFLU.PROG, IFLULA.PROG, CELLT,
766 : & LOGAN)
767 : ELSEIF(INDMET .EQ. 3)THEN
768 : C
769 : C********** AUSMP-up
770 : C
771 : V_INF=MAX(MPUIN1.VPOCHA(NLCEG,1),
772 : & MPUIN1.VPOCHA(NLCED,1))
772 : CALL FAUSRE(IDIM, NESP, NLHS, NORD,
773 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
774 : & TMAX,V_INF,
775 : & CNNA, CTNA, CVNA,
776 : & IWST_L.PROG,
777 : & IWST_R.PROG,
778 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
779 : & ICVRSS.PROG, ICVWOR.PROG,
780 : & IFLU.PROG, IFLULA.PROG, CELLT,
781 : & LOGAN)
782 : ENDIF
783 :
784 : C write(*,*) 'Riemann problem 1-1'
785 : C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
786 : C
787 : C subroutine fcolre(ndim, nesp, nLHS, nordpo,
788 : C & reacoe, aprop, Runiv,
789 : C & Tmaxcv,
790 : C & wvect_l,
791 : C & wvect_r,
792 : C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
793 : C & flu, flu_lag, cellt,
794 : C & logan)
795 : C
796 : C
797 : IF (LOGAN) THEN
798 : WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
799 : WRITE(IOIMP,*) 'RIEMANN PROBLEM 1-1'
800 : WRITE(IOIMP,*) 'FACE ', NGCF, ' S11 ', S1_11
801 : WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
802 : DO I2 = 1, 2 + IDIM
803 : WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
804 : $ IWST_R.PROG(I2)
805 : ENDDO
806 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
807 : GOTO 9999
808 : ENDIF
809 : C
810 : CELLTM = MAX(CELLT, CELLTM)
811 : C
812 : C************* Ecriture du residu
813 : C
814 : C FLUX(2) = RO Un RO Un
815 : C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
816 : C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
817 : C FLUX(2 + IDIM) = RO Un Et RO Un Et
818 : C
819 : SURF = SURF0 * S1_11
820 : FLUCEL = IFLU.PROG(1) * SURF
821 : C write(*,*) 'flux', iflu.prog(1)
822 : C write(*,*) 'nlceg, nlced', nlceg, nlced
823 : C write(*,*) 'surf', surf
824 : C
825 : C Rho un
826 : C
827 : MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
828 : & (FLUCEL / VOLG)
829 : IF (NLCED .NE. NLCEG)
830 : & MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
831 : & (FLUCEL / VOLD)
832 : C
833 : C Rho (Un Ux + P CNX)
834 : C NB : In 2D (IFLU.PROG(4)*CVX) = RO Un Et * 0 = 0
835 : C In 3D (IFLU.PROG(4)*CVX) = RO Un uv * CVX
836 : C
837 : FLUCEL = ( (IFLU.PROG(2)*CNX) + (IFLU.PROG(3)*CTX) +
838 : & (IFLU.PROG(4)*CVX)) * SURF
839 : MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
840 : & (FLUCEL / VOLG)
841 : IF (NLCED .NE. NLCEG)
842 : & MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
843 : & (FLUCEL / VOLD)
844 : C
845 : C Rho (Un Uy + P CNY)
846 : C
847 : FLUCEL = ( (IFLU.PROG(2)*CNY) + (IFLU.PROG(3)*CTY)
848 : & + (IFLU.PROG(4)*CVY) ) * SURF
849 : MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
850 : & (FLUCEL / VOLG)
851 : IF (NLCED .NE. NLCEG)
852 : & MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
853 : & (FLUCEL / VOLD)
854 : C
855 : C Rho (Un Uz + P CNZ)
856 : C
857 : IF (IDIM .EQ. 3) THEN
858 : FLUCEL = ( (IFLU.PROG(2)*CNZ) + (IFLU.PROG(3)*CTZ)
859 : & + (IFLU.PROG(4)*CVZ) ) * SURF
860 : MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
861 : & (FLUCEL / VOLG)
862 : IF (NLCED .NE. NLCEG)
863 : & MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
864 : & (FLUCEL / VOLD)
865 : ENDIF
866 : C
867 : C Rho Un et
868 : C
869 : FLUCEL = IFLU.PROG(2 + IDIM) * SURF
870 : MPORES.VPOCHA(NLCEG,3 + IDIM) =
871 : & MPORES.VPOCHA(NLCEG, 3 + IDIM) -
872 : & (FLUCEL / VOLG)
873 : IF (NLCED .NE. NLCEG)
874 : & MPORES.VPOCHA(NLCED, 3 + IDIM) =
875 : & MPORES.VPOCHA(NLCED, 3 + IDIM) +
876 : & (FLUCEL / VOLD)
877 : C
878 : ENDIF
879 : C
880 : C********** Riemann problem 2-2
881 : C
882 : IF (S2_22 .GE. EPS) THEN
883 : IWST_L.PROG(1) = ROG2
884 : IWST_L.PROG(2) = UNG2
885 : C IWST_L.PROG(2) = 0.0D0
886 : IWST_L.PROG(3) = UTG2
887 : IWST_L.PROG(IDIM + 2) = PG2
888 : IWST_L.PROG(IDIM + 3) = 1.0D0
889 : C
890 : IWST_R.PROG(1) = ROD2
891 : IWST_R.PROG(2) = UND2
892 : C IWST_R.PROG(2) = 0.0D0
893 : IWST_R.PROG(3) = UTD2
894 : IWST_R.PROG(IDIM + 2) = PD2
895 : IWST_R.PROG(IDIM + 3) = 1.0D0
896 : IF (IDIM .EQ. 3) THEN
897 : IWST_L.PROG(4) = UVG2
898 : IWST_R.PROG(4) = UVD2
899 : ENDIF
900 : C
901 : C do i2 = 1 , 4 + NESP + IDIM, 1
902 : C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
903 : C enddo
904 : C
905 : IF (INDMET .EQ. 1) THEN
905 : CALL FVLHRE(IDIM, NESP, NLHS, NORD,
906 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
907 : & TMAX,
908 : & CNNA, CTNA, CVNA,
909 : & IWST_L.PROG,
910 : & IWST_R.PROG,
911 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
912 : & ICVRSS.PROG, ICVWOR.PROG,
913 : & IFLU.PROG, IFLULA.PROG, CELLT,
914 : & LOGAN)
915 : ELSEIF(INDMET .EQ. 2)THEN
915 : CALL FBECRE(IDIM, NESP, NLHS, NORD,
916 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
917 : & TMAX,
918 : & CNNA, CTNA, CVNA,
919 : & IWST_L.PROG,
920 : & IWST_R.PROG,
921 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
922 : & ICVRSS.PROG, ICVWOR.PROG,
923 : & IFLU.PROG, IFLULA.PROG, CELLT,
924 : & LOGAN)
925 : ELSEIF(INDMET .EQ. 3)THEN
926 : C
927 : C********** AUSMP-up
928 : C
929 : V_INF=MAX(MPUIN2.VPOCHA(NLCEG,1),
930 : & MPUIN2.VPOCHA(NLCED,1))
930 : CALL FAUSRE(IDIM, NESP, NLHS, NORD,
931 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
932 : & TMAX, V_INF,
933 : & CNNA, CTNA, CVNA,
934 : & IWST_L.PROG,
935 : & IWST_R.PROG,
936 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
937 : & ICVRSS.PROG, ICVWOR.PROG,
938 : & IFLU.PROG, IFLULA.PROG, CELLT,
939 : & LOGAN)
940 : ENDIF
941 : C
942 : C write(*,*) 'Riemann problem 2-2'
943 : C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
944 : C
945 : C subroutine fcolre(ndim, nesp, nLHS, nordpo,
946 : C & reacoe, aprop, Runiv,
947 : C & Tmaxcv,
948 : C & wvect_l,
949 : C & wvect_r,
950 : C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
951 : C & flu, flu_lag, cellt,
952 : C & logan)
953 : C
954 : C
955 : IF (LOGAN) THEN
956 : WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
957 : WRITE(IOIMP,*) 'RIEMANN PROBLEM 2-2'
958 : WRITE(IOIMP,*) 'FACE ', NGCF, ' S22 ', S2_22
959 : WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
960 : DO I2 = 1, 2 + IDIM
961 : WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
962 : $ IWST_R.PROG(I2)
963 : ENDDO
964 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
965 : GOTO 9999
966 : ENDIF
967 : C
968 : CELLTM = MAX(CELLT, CELLTM)
969 : C
970 : C
971 : C************* Ecriture du residu
972 : C
973 : C FLUX(7) = RO Un RO Un
974 : C FLUX(8) = RO Un Un + P -> RO Un Ux + P CNX
975 : C FLUX(9) = RO Un Ut -> RO Un Uy + P CNY
976 : C FLUX(10)= RO Un Et RO Un Et
977 : C
978 : SURF = SURF0 * S2_22
979 : C
980 : C Rho un
981 : C
982 : FLUCEL = IFLU.PROG(1) * SURF
983 : MPORES.VPOCHA(NLCEG, (IDIM + 5)) =
984 : & MPORES.VPOCHA(NLCEG, (IDIM + 5)) -
985 : & (FLUCEL / VOLG)
986 : IF (NLCED .NE. NLCEG)
987 : & MPORES.VPOCHA(NLCED, (IDIM + 5)) =
988 : & MPORES.VPOCHA(NLCED, (IDIM + 5)) +
989 : & (FLUCEL / VOLD)
990 : C
991 : C Rho un ux + P cnx
992 : C
993 : FLUCEL = ((IFLU.PROG(2)*CNX) + (IFLU.PROG(3)*CTX) +
994 : & (IFLU.PROG(4)*CVX)) * SURF
995 : MPORES.VPOCHA(NLCEG, (IDIM + 6)) =
996 : & MPORES.VPOCHA(NLCEG, (IDIM + 6)) -
997 : & (FLUCEL / VOLG)
998 : IF (NLCED .NE. NLCEG)
999 : & MPORES.VPOCHA(NLCED, (IDIM + 6)) =
1000 : & MPORES.VPOCHA(NLCED, (IDIM + 6)) +
1001 : & (FLUCEL / VOLD)
1002 : C
1003 : C Rho un uy + P cny
1004 : C
1005 : FLUCEL = ((IFLU.PROG(2)*CNY) + (IFLU.PROG(3)*CTY) +
1006 : & (IFLU.PROG(4)*CVY)) * SURF
1007 : MPORES.VPOCHA(NLCEG, (IDIM + 7)) =
1008 : & MPORES.VPOCHA(NLCEG, (IDIM + 7)) -
1009 : & (FLUCEL / VOLG)
1010 : IF (NLCED .NE. NLCEG)
1011 : & MPORES.VPOCHA(NLCED, (IDIM + 7)) =
1012 : & MPORES.VPOCHA(NLCED, (IDIM + 7)) +
1013 : & (FLUCEL / VOLD)
1014 : C
1015 : C Rho un uz + P cnz
1016 : C
1017 : IF (IDIM .EQ. 3) THEN
1018 : FLUCEL = ((IFLU.PROG(2)*CNZ) + (IFLU.PROG(3)*CTZ) +
1019 : & (IFLU.PROG(4)*CVZ)) * SURF
1020 : MPORES.VPOCHA(NLCEG, (IDIM + 8)) =
1021 : & MPORES.VPOCHA(NLCEG, (IDIM + 8)) -
1022 : & (FLUCEL / VOLG)
1023 : IF (NLCED .NE. NLCEG)
1024 : & MPORES.VPOCHA(NLCED, (IDIM + 8)) =
1025 : & MPORES.VPOCHA(NLCED, (IDIM + 8)) +
1026 : & (FLUCEL / VOLD)
1027 : ENDIF
1028 : C
1029 : C Rho un et
1030 : C
1031 : FLUCEL = IFLU.PROG(IDIM + 2) * SURF
1032 : MPORES.VPOCHA(NLCEG, ((2 * IDIM) + 6)) =
1033 : & MPORES.VPOCHA(NLCEG, ((2 * IDIM) + 6)) -
1034 : & (FLUCEL / VOLG)
1035 : IF (NLCED .NE. NLCEG)
1036 : & MPORES.VPOCHA(NLCED, ((2 * IDIM) + 6)) =
1037 : & MPORES.VPOCHA(NLCED, ((2 * IDIM) + 6)) +
1038 : & (FLUCEL / VOLD)
1039 : ENDIF
1040 : C
1041 : C********** Riemann problem 1-2
1042 : C From a numerical point of view, S1_12 and S2_12 are identical.
1043 : C The volume fraction of 1 is larger in l than in r.
1044 : C
1045 : C rho, u, p for the phase 1 and 2.
1046 : C NB in case of wall-boundary conditions, we cannot be here
1047 : C since ALPD1 = ALPG1, ALPD2 = ALPG2.
1048 : C
1049 : IF ((S1_12 .gt. EPS) .or. (S2_12 .gt. EPS)) THEN
1050 : C
1051 : IWST_L.PROG(1) = ROG1
1052 : IWST_L.PROG(2) = UNG1
1053 : C IWST_L.PROG(2) = 0.0D0
1054 : IWST_L.PROG(3) = UTG1
1055 : IWST_L.PROG(IDIM + 2) = PG1
1056 : IWST_L.PROG(IDIM + 3) = 0.0D0
1057 : C
1058 : IWST_R.PROG(1) = ROD2
1059 : IWST_R.PROG(2) = UND2
1060 : C IWST_R.PROG(2) = 0.0D0
1061 : IWST_R.PROG(3) = UTD2
1062 : IWST_R.PROG(IDIM + 2) = PD2
1063 : IWST_R.PROG(IDIM + 3) = 1.0D0
1064 : IF (IDIM .EQ. 3) THEN
1065 : IWST_L.PROG(4) = UVG1
1066 : IWST_R.PROG(4) = UVD2
1067 : ENDIF
1068 : C
1069 : C do i2 = 1 , 4 + NESP + IDIM, 1
1070 : C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
1071 : C enddo
1072 : C
1072 : CALL FBECRE(IDIM, NESP, NLHS, NORD,
1073 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
1074 : & TMAX,
1075 : & CNNA, CTNA, CVNA,
1076 : & IWST_L.PROG,
1077 : & IWST_R.PROG,
1078 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
1079 : & ICVRSS.PROG, ICVWOR.PROG,
1080 : & IFLU.PROG, IFLULA.PROG, CELLT,
1081 : & LOGAN)
1082 : C
1083 : C write(*,*) 'Riemann problem 1-2'
1084 : C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
1085 : C
1086 : C subroutine fcolre(ndim, nesp, nLHS, nordpo,
1087 : C & reacoe, aprop, Runiv,
1088 : C & Tmaxcv,
1089 : C & wvect_l,
1090 : C & wvect_r,
1091 : C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
1092 : C & flu, flu_lag, cellt,
1093 : C & logan)
1094 : C
1095 : IF (LOGAN) THEN
1096 : WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
1097 : WRITE(IOIMP,*) 'RIEMANN PROBLEM 1-2'
1098 : WRITE(IOIMP,*) 'FACE ', NGCF, ' S12 ', S1_12, S2_12
1099 : WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
1100 : DO I2 = 1, 2 + IDIM
1101 : WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
1102 : $ IWST_R.PROG(I2)
1103 : ENDDO
1104 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
1105 : GOTO 9999
1106 : ENDIF
1107 : C
1108 : CELLTM = MAX(CELLT, CELLTM)
1109 : C
1110 : C************* Ecriture du residu
1111 : C
1112 : C FLUX(7) = RO Un RO Un
1113 : C FLUX(8) = RO Un Un + P -> RO Un Ux + P CNX
1114 : C FLUX(9) = RO Un Ut -> RO Un Uy + P CNY
1115 : C FLUX(10)= RO Un Et RO Un Et
1116 : C
1117 : C FLULA(3+IDIM) = - D_resh, with D_resh is the reactive shock
1118 : C speed
1119 : C
1120 : C We update resu (we are dealing with a reactive
1121 : C Riemann problem)
1122 : C
1123 : IF (IFLULA.PROG(IDIM + 3) .LT. 0.0D0) THEN
1124 : C
1125 : C Note that IFLULA.PROG(IDIM + 3) = - D_resh
1126 : C
1127 : C t
1128 : C |F / Flag
1129 : C ---> --->
1130 : C | /
1131 : C 1 | / 2
1132 : C | /
1133 : C |/
1134 : C ------------------------------> x
1135 : C ifac ifac+1
1136 : C
1137 : C D_resh > 0
1138 : C IFLULA.PROG(5) < 0
1139 : C
1140 : C Definition of SCON (conservative surface) and SG and SD
1141 : C (lagrangian surfaces)
1142 : C
1143 : c NB SCON > 0 if, for a positive conservative flux, gives negative
1144 : C contribution in ifac and positive contribution
1145 : C in ifac+1
1146 : C If follows that SCON \ge 0
1147 : C SG and SD > 0 if, for a positive lagrangian flux, gives
1148 : C negative contribution in ifac and negative
1149 : C contributuion in ifac+1
1150 : C
1151 :
1152 : C In first order case :
1153 : C PHASE 1, conservative contribution in ifac (-F) and ifac+1 (+F)
1154 : C non-conservative contribution in ifac+1 (-Flag)
1155 : C negative contribution => SG and SD > 0
1156 : C
1157 : SCON1 = S1_12E
1158 : SD1 = S1_12E + S1_12R
1159 : SG1 = S1_12L
1160 : C In first order case :
1161 : C PHASE 2, no conservative contribution
1162 : C non-conservative contribution in ifac+1 (Flag)
1163 : C positive contribution => SG and SD < 0
1164 : C
1165 : SCON2 = 0.0D0
1166 : SD2 = -1.0D0 * (S2_12E + S2_12R)
1167 : SG2 = -1.0D0 * S2_12L
1168 : ELSE
1169 : C
1170 : C t
1171 : C \Flag | F
1172 : C ---> --->
1173 : C \ |
1174 : C 1 \ | 2
1175 : C \ |
1176 : C \|
1177 : C ------------------------------> x
1178 : C ifac ifac+1
1179 : C
1180 : C D_resh < 0
1181 : C IFLULA.PROG(IDIM + 3) > 0
1182 : C
1183 : C In first order case :
1184 : C PHASE 1, no conservative contribution in ifac and ifac+1
1185 : C non-conservative contribution in ifac (-Flag)
1186 : C negative contribution => SG and SD > 0
1187 : C
1188 : SCON1 = 0.0D0
1189 : SG1 = S1_12E + S1_12L
1190 : SD1 = S1_12R
1191 : C In first order case :
1192 : C PHASE 2, conservative contribution in ifac and ifac+1
1193 : C non-conservative contribution in ifac (Flag)
1194 : C positive contribution => SG and SD < 0
1195 : C
1196 : SCON2 = S2_12E
1197 : SG2 = -1.0D0 * (S2_12L + S2_12E)
1198 : SD2 = -1.0D0 * S2_12R
1199 : ENDIF
1200 : SURF = SURF0
1200 : CALL RDEMUP(IDIM, 0,
1201 : & SCON1, SG1, SD1,
1202 : & SCON2, SG2, SD2,
1203 : & SURF,VOLG,VOLD,
1204 : & CNX,CNY,CNZ,CTX,CTY,CTZ,CVX,CVY,CVZ,
1205 : & IFLU.PROG,IFLULA.PROG,RESG,RESD)
1206 : C write(*,*) 'RESD'
1207 : DO ICOMP = 1, ((2 * IDIM) + 6)
1208 : MPORES.VPOCHA(NLCED,ICOMP) =
1209 : & MPORES.VPOCHA(NLCED,ICOMP) + RESD(ICOMP)
1210 : C write(*,*) ICOMP, RESD(ICOMP)
1211 : ENDDO
1212 : C write(*,*) 'RESG'
1213 : DO ICOMP = 1, ((2 * IDIM) + 6)
1214 : MPORES.VPOCHA(NLCEG,ICOMP) =
1215 : & MPORES.VPOCHA(NLCEG,ICOMP) + RESG(ICOMP)
1216 : C write(*,*) ICOMP, RESG(ICOMP)
1217 : ENDDO
1218 : ENDIF
1219 : C
1220 : C********** Riemann problem 2-1
1221 : C From a numerical point of view, S2_21 and S2_21 are identical.
1222 : C The volume fraction of 2 is larger in l than in r.
1223 : C
1224 : C rho, u, p for the phase 1 and 2.
1225 : C
1226 : C NB in case of wall-boundary conditions, we cannot be here
1227 : C since ALPD1 = ALPG1, ALPD2 = ALPG2.
1228 : C
1229 : IF ((S1_21 .gt. EPS) .or. (S2_21 .gt. EPS)) THEN
1230 : C
1231 : C
1232 : IWST_L.PROG(1) = ROG2
1233 : IWST_L.PROG(2) = UNG2
1234 : CC IWST_L.PROG(2) = 0.0D0
1235 : IWST_L.PROG(3) = UTG2
1236 : IWST_L.PROG(IDIM + 2) = PG2
1237 : IWST_L.PROG(IDIM + 3) = 1.0D0
1238 : CC
1239 : IWST_R.PROG(1) = ROD1
1240 : IWST_R.PROG(2) = UND1
1241 : C IWST_R.PROG(2) = 0.0D0
1242 : IWST_R.PROG(3) = UTD1
1243 : IWST_R.PROG(IDIM + 2) = PD1
1244 : IWST_R.PROG(IDIM + 3) = 0.0D0
1245 : C
1246 : IF (IDIM .EQ. 3) THEN
1247 : IWST_L.PROG(4) = UVG2
1248 : IWST_R.PROG(4) = UVD1
1249 : ENDIF
1250 : C
1251 : CC
1252 : C do i2 = 1 , 4 + NESP + IDIM, 1
1253 : C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
1254 : C enddo
1255 : C
1255 : CALL FBECRE(IDIM, NESP, NLHS, NORD,
1256 : & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
1257 : & TMAX,
1258 : & CNNA, CTNA, CVNA,
1259 : & IWST_L.PROG,
1260 : & IWST_R.PROG,
1261 : & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
1262 : & ICVRSS.PROG, ICVWOR.PROG,
1263 : & IFLU.PROG, IFLULA.PROG, CELLT,
1264 : & LOGAN)
1265 : C write(*,*) 'Riemann problem 1-2'
1266 : C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
1267 : C
1268 : C subroutine fcolre(ndim, nesp, nLHS, nordpo,
1269 : C & reacoe, aprop, Runiv,
1270 : C & Tmaxcv,
1271 : C & wvect_l,
1272 : C & wvect_r,
1273 : C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
1274 : C & flu, flu_lag, cellt,
1275 : C & logan)
1276 : C
1277 : IF (LOGAN) THEN
1278 : WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
1279 : WRITE(IOIMP,*) 'RIEMANN PROBLEM 2-1'
1280 : WRITE(IOIMP,*) 'FACE ', NGCF, ' S21 ', S1_21, S2_21
1281 : WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
1282 : DO I2 = 1, 2 + IDIM
1283 : WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
1284 : $ IWST_R.PROG(I2)
1285 : ENDDO
1286 : WRITE(IOIMP,*) 'ANOMALY DETECTED'
1287 : GOTO 9999
1288 : ENDIF
1289 : C
1290 : CELLTM = MAX(CELLT, CELLTM)
1291 : C
1292 : C************* Ecriture du residu
1293 : C
1294 : C FLUX(7) = RO Un RO Un
1295 : C FLUX(8) = RO Un Un + P -> RO Un Ux + P CNX
1296 : C FLUX(9) = RO Un Ut -> RO Un Uy + P CNY
1297 : C FLUX(10)= RO Un Et RO Un Et
1298 : C
1299 : C FLULA(3+IDIM) = - D_resh, with D_resh is the reactive shock
1300 : C speed
1301 : C
1302 : C
1303 : C We update resu (we are dealing with a reactive
1304 : C Riemann problem)
1305 : C
1306 : IF (IFLULA.PROG(IDIM + 3) .LT. 0.0D0) THEN
1307 : C
1308 : C Note that IFLULA.PROG(IDIM + 3) = - D_resh
1309 : C
1310 : C t
1311 : C |F / Flag
1312 : C ---> --->
1313 : C | /
1314 : C 2 | / 1
1315 : C | /
1316 : C |/
1317 : C ------------------------------> x
1318 : C ifac ifac+1
1319 : C
1320 : C D_resh > 0
1321 : C IFLULA.PROG(IDIM + 3) < 0
1322 : C
1323 : C In first order case :
1324 : C PHASE 2, conservative contribution in ifac (-F) and ifac+1 (+F)
1325 : C non-conservative contribution in ifac+1 (-Flag)
1326 : C negative contribution => SG and SD > 0
1327 : C
1328 : SCON2 = S2_21E
1329 : SG2 = S2_21L
1330 : SD2 = S2_21E + S2_21R
1331 : C
1332 : C PHASE 1, non-conservative contribution in ifac+1 (+Flag)
1333 : C positive contribution => SG and SD < 0
1334 : C
1335 : SCON1 = 0.0D0
1336 : SG1 = -1.0D0 * S1_21L
1337 : SD1 = -1.0D0 * (S1_21E + S1_21R)
1338 : C
1339 : ELSE
1340 : C
1341 : C t
1342 : C \Flag | F
1343 : C ---> --->
1344 : C \ |
1345 : C 2 \ | 1
1346 : C \ |
1347 : C \|
1348 : C ------------------------------> x
1349 : C ifac ifac+1
1350 : C
1351 : C D_resh < 0
1352 : C IFLULA.PROG(IDIM + 3) > 0
1353 : C
1354 : C In first order case :
1355 : C PHASE 2, no conservative contribution
1356 : C non-conservative contribution in ifac (-Flag)
1357 : C negative contribution => SG and SD > 0
1358 : C
1359 : SCON2 = 0.0D0
1360 : SG2 = S2_21E + S2_21L
1361 : SD2 = S2_21R
1362 : C
1363 : C In first order case :
1364 : C PHASE 1, conservative contribution in ifac (-F) and
1365 : C ifac+1 (+F)
1366 : C non-conservative contribution in ifac (+Flag)
1367 : C positive contribution => SG and SD > 0
1368 : C
1369 : C ALPHA1
1370 : C
1371 : SCON1 = S1_21E
1372 : SG1 = -1.0D0 * (S1_21E + S1_21L)
1373 : SD1 = -1.0D0 * S1_21R
1374 : ENDIF
1375 : SURF = SURF0
1375 : CALL RDEMUP(IDIM, 0,
1376 : & SCON1, SG1, SD1,
1377 : & SCON2, SG2, SD2,
1378 : & SURF,VOLG,VOLD,
1379 : & CNX,CNY,CNZ,CTX,CTY,CTZ,CVX,CVY,CVZ,
1380 : & IFLU.PROG,IFLULA.PROG,RESG,RESD)
1381 : C write(*,*) 'RESD'
1382 : DO ICOMP = 1, ((2 * IDIM) + 6)
1383 : MPORES.VPOCHA(NLCED,ICOMP) =
1384 : & MPORES.VPOCHA(NLCED,ICOMP) + RESD(ICOMP)
1385 : C write(*,*) ICOMP, RESD(ICOMP)
1386 : ENDDO
1387 : C write(*,*) 'RESG'
1388 : DO ICOMP = 1, ((2 * IDIM) + 6)
1389 : MPORES.VPOCHA(NLCEG,ICOMP) =
1390 : & MPORES.VPOCHA(NLCEG,ICOMP) + RESG(ICOMP)
1391 : C write(*,*) ICOMP, RESG(ICOMP)
1392 : ENDDO
1393 : ENDIF
1394 : C
1395 : C********** Calcul du pas du temps (CFL)
1396 : C
1397 : C********* a) etat a l'interface
1398 : C
1399 : DIAMG = MPOVDI.VPOCHA(NLCEG,1)
1400 : DIAMD = MPOVDI.VPOCHA(NLCED,1)
1401 : DIAM = MIN(DIAMG,DIAMD)
1402 : CELL = CELLTM/DIAM
1403 : IF(CELL .GT. UNSDT)THEN
1404 : UNSDT = CELL
1405 : ENDIF
1406 : ENDIF
1407 : C
1408 : C**** Fin boucle sur FACEL
1409 : C
1410 : ENDDO
1411 : C
1412 : C**** Pas du temps (condition de non interaction en 1D)
1413 : C
1414 : DT = 0.5D0 / UNSDT
1415 : C
1416 : C**** Desactivation des segments
1417 : C
1418 : SEGSUP MLELIM
1419 : SEGSUP MLENT1
1420 : SEGSUP MLENT2
1421 : C
1422 : SEGDES MPOVSU
1423 : SEGDES MPOVDI
1424 : SEGDES MPOVVO
1425 : SEGDES MPORES
1426 : SEGDES IPT2
1427 : C
1428 : C Warning. MLRMFR = PGAS . 'MASSFRA'
1429 : C MLRCHE = PGAS . 'CHEMCOEF'
1430 : C
1431 : SEGDES MLRMFR
1432 : SEGDES MLRCHE
1433 : SEGSUP IY
1434 : C
1435 : SEGSUP IWWORK
1436 : SEGSUP IWST_R
1437 : SEGSUP IWST_L
1438 : SEGSUP ICVL
1439 : SEGSUP ICVR
1440 : SEGSUP ICVRSS
1441 : SEGSUP ICVWOR
1442 : SEGSUP IFLU
1443 : SEGSUP IFLULA
1444 : C
1445 : SEGDES MPALP1
1446 : SEGDES MPALP2
1447 : SEGDES MELAL1
1448 : SEGDES MELRO1
1449 : SEGDES MELP1
1450 : SEGDES MELAL2
1451 : SEGDES MELRO2
1452 : SEGDES MELP2
1453 : SEGDES MELUN2
1454 : SEGDES MELUT2
1455 : IF (IDIM .EQ. 2) THEN
1456 : SEGDES MELUN1
1457 : SEGDES MELUT1
1458 : SEGDES MELVNX
1459 : SEGDES MELVNY
1460 : SEGDES MELT1X
1461 : SEGDES MELT1Y
1462 : SEGDES MELUN2
1463 : SEGDES MELUT2
1464 : ELSE
1465 : SEGDES MELVNX
1466 : SEGDES MELVNY
1467 : SEGDES MELVNZ
1468 : SEGDES MELT1X
1469 : SEGDES MELT1Y
1470 : SEGDES MELT1Z
1471 : SEGDES MELT2X
1472 : SEGDES MELT2Y
1473 : SEGDES MELT2Z
1474 : SEGDES MELUN1
1475 : SEGDES MELUT1
1476 : SEGDES MELUV1
1477 : SEGDES MELUN2
1478 : SEGDES MELUT2
1479 : SEGDES MELUV2
1480 : ENDIF
1481 : C
1482 : SEGDES MPOALP
1483 : IF (INDK0 .EQ. 2) SEGDES MPOK0
1484 : C
1485 : IF(MPUIN2 .NE. 0)THEN
1486 : SEGDES MPUIN2
1487 : SEGDES MPUIN1
1488 : ENDIF
1489 : C
1490 : 9999 CONTINUE
1491 : RETURN
1492 : END
1493 : CC
1494 :
1495 :
1496 :
1497 :
© Cast3M 2003 - All rights reserved.
Disclaimer