1 : C PRE221 SOURCE BECC 10/10/07 21:15:00 6774 2 : SUBROUTINE PRE221(LOGTEM, 3 : & ICEN,IFACE,IFACEL,INORM, 4 : & IROC, IGRROC, IALROC, 5 : & IVITC, IGRVC, IALVC, 6 : & IPC ,IGRPC, IALPC, 7 : & MLRECV, MLRECP, MLMOT1, 8 : & IYC ,IGRYC, IALYC, 9 : & IGAMC, 10 : & DELTAT, 11 : & IROF,IVITF,IPF,IGAMF,IYF, 12 : & LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2) 13 : C************************************************************************ 14 : C 15 : C PROJET : CASTEM 2000 16 : C 17 : C NOM : PRE221 18 : C 19 : C DESCRIPTION : Voir PRE22 20 : C 21 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) 22 : C 23 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF 24 : C 25 : C************************************************************************ 26 : C 27 : C 28 : C APPELES (Outils) : KRIPAD, LICHT 29 : C 30 : C APPELES (Calcul) : AUCUN 31 : C 32 : C 33 : C************************************************************************ 34 : C 35 : C ENTREES 36 : C 37 : C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps 38 : C sinon 1er ordre en temps 39 : C 40 : C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE 41 : C 42 : C ICEN : MELEME de 'POI1' SPG des CENTRES 43 : C 44 : C IFACE : MELEME de 'POI1' SPG des FACES 45 : C 46 : C IFACEL : MELEME de 'SEG3' avec 47 : C CENTRE d'Elt "gauche" 48 : C CENTRE de Face 49 : C CENTRE d'Elt "droite" 50 : C 51 : C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2) 52 : C 53 : C INORM : CHPOINT des cosinus directeurs de normales aux faces 54 : C 55 : C 2) Autres pointeurs 56 : C 57 : C IROC : CHPOINT "CENTRE" contenant la masse volumique RHO 58 : C 59 : C IGRROC : CHPOINT "CENTRE" contenant le gradient de la 60 : C masse volumique RHO (2 composantes) 61 : C 62 : C IALROC : CHPOINT "CENTRE" contenant le limiteur du gradient 63 : C de la masse volumique 64 : C 65 : C IVITC : CHPOINT "CENTRE" contenant la vitesse UX, UY ; 66 : C 67 : C IGRVC : CHPOINT "CENTRE" contenant le gradient de la 68 : C vitesse (4 composantes) 69 : C 70 : C IALVC : CHPOINT "CENTRE" contenant le limiteur du gradient 71 : C de la vitesse (2 composantes) 72 : C 73 : C IPC : CHPOINT "CENTRE" contenat la pression P; 74 : C 75 : C IGRPC : CHPOINT "CENTRE" contenant le gradient de la 76 : C pression (2 composantes) 77 : C 78 : C IALPC : CHPOINT "CENTRE" contenant le limiteur du gradient 79 : C de la pression 80 : C 81 : C MLRECV : pointeur de la LISTREEL des cv 82 : C 83 : C MLRECP : pointeur de la LISTREEL des cp 84 : C 85 : C MLMOT1 : pointeur des noms des gaz qui apparessent dans les 86 : C equations d'Euler 87 : C 88 : C IYC : CHPOINT "CENTRE" contenat les fractions massiques ; 89 : C 90 : C IGRYC : CHPOINT "CENTRE" contenant les gradient des fr.mass.; 91 : C 92 : C IALPC : CHPOINT "CENTRE" contenant les limiteurs des gradients 93 : C des fr.mass.; 94 : C 95 : C IGAMC : CHPOINT "CENTRE" contenat le "Gamma" du gaz 96 : C 97 : C 3) 98 : C 99 : C DELTAT : REAL*8, encrement en temps pour calculer la prediction 100 : C 101 : C 102 : C SORTIES 103 : C 104 : C 105 : C IROF : MCHAML defini sur le MELEME de pointeur IFACEL, 106 : C contenant la masse volumique RHO 107 : C 108 : C IVITF : MCHAML defini sur le MELEME de pointeur IFACEL, 109 : C contenant la vitesse UN, UT dans le repaire local 110 : C (n,t) et defini sur le MELEME de pointeur IFACE, 111 : C contenant les cosinus directeurs du repere local 112 : C 113 : C IPF : MCHAML defini sur le MELEME de pointeur IFACEL, 114 : C contenant la pression P 115 : C 116 : C IYF : MCHAML defini sur le MELEME de pointeur IFACEL, 117 : C contenant les fractions massiques; 118 : C 119 : C IGAMF : MCHAML defini sur le MELEME de pointeur IFACEL, 120 : C contenant le "gamma" du gaz 121 : C 122 : C LOGAN : anomalie detectee 123 : C 124 : C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité 125 : C negative a été detectée -> en interactif le 126 : C programme s'arrete en GIBIANE 127 : C (erreur stocké en MESERR et VALER) 128 : C 129 : C LOGBOR : (LOGICAL): si .TRUE. un gamma a ete detecte 130 : C dehor 1 et 3 (sa valeur stockée en MESERR et VALER; 131 : C en VAL1 et en VAL2 on stocke 1.0 et 3.0) 132 : C 133 : C MESERR 134 : C VALER 135 : C VAL1, 136 : C VAL2 : pour les messages d'erreur 137 : C 138 : C************************************************************************ 139 : C 140 : C HISTORIQUE (Anomalies et modifications éventuelles) 141 : C 142 : C HISTORIQUE : Créée le 10.7.98. 143 : C 144 : C************************************************************************ 145 : C 146 : C 147 : C ATTENTION: Cet programme marche si le MAILLAGE est convex; 148 : C si non il faut changer l'algoritme de calcul de 149 : C l'orientation des normales aux faces. 150 : C 151 : C La positivité n'est pas controlle parce que c'est déjà fait 152 : C dans l'operateur PRIM 153 : C 154 : C 155 : C************************************************************************ 156 : C 157 : C**** Variables de COOPTIO 158 : C 159 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 160 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 161 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 162 : C & ,IECHO, IIMPI, IOSPI 163 : C & ,IDIM 164 : CC & ,MCOORD 165 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 166 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 167 : C & ,NORINC,NORVAL,NORIND,NORVAD 168 : C & ,NUCROU, IPSAUV 169 : C 170 : C**** Les variables 171 : C 172 : IMPLICIT INTEGER(I-N) 173 : INTEGER ICEN, IFACE, IFACEL, INORM,IROC, IGRROC, IALROC 174 : & , IVITC, IGRVC, IALVC 175 : & , IPC ,IGRPC, IALPC 176 : & , IYC, IGRYC, IALYC 177 : & , IGAMC 178 : & , IROF, IVITF, IPF, IGAMF, IYF 179 : & , IGEOM, NFAC, NCEN 180 : & , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE 181 : & , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1 182 : & , IDIMP1, INDCEL, I1, NESP 183 : REAL*8 VALER, VAL1, VAL2, XG, YG, XC, YC, XD, YD, DELTAT 184 : & ,DXG, DYG, DXD, DYD 185 : & , CNX, CNY, CTX, CTY, ORIENT 186 : & , ROG, PG, GAMG, UXG, UYG, UNG, UTG 187 : & , ROD, PD, GAMD, UXD, UYD, UND, UTD 188 : & , VALCEL, DCEL, ALCEL 189 : & , DROX, DROY, DUXX, DUXY, DUYX, DUYY, DPX, DPY 190 : & , DRO, DUX, DUY, DP 191 : & , ALPHA, CPCEL, CVCEL, DYMAS, SUMY 192 : CHARACTER*(40) MESERR 193 : CHARACTER*(8) TYPE 194 : CHARACTER*(4) CARCEL 195 : LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGI2 196 : C 197 : C 198 : C**** Les Includes 199 : C 200 : -INC SMCOORD 201 : -INC CCOPTIO 202 : -INC SMCHPOI 203 : POINTEUR MPROC.MPOVAL, MPGRR.MPOVAL, 204 : & MPVITC.MPOVAL, MPGRV.MPOVAL, 205 : & MPPC.MPOVAL, MPGRP.MPOVAL, 206 : & MPYC.MPOVAL, MPGRY.MPOVAL, 207 : & MPGAMC.MPOVAL, MPNORM.MPOVAL, 208 : & MPROP.MPOVAL, MPPP.MPOVAL, MPVITP.MPOVAL, 209 : & MPYP.MPOVAL 210 : -INC SMCHAML 211 : POINTEUR MCHAMY.MCHAML 212 : POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, 213 : & MELT1X.MELVAL, MELT1Y.MELVAL 214 : POINTEUR MELVUN.MELVAL, MELVUT.MELVAL 215 : POINTEUR MELRO.MELVAL, MELP.MELVAL, 216 : & MELGAM.MELVAL 217 : POINTEUR MELALR.MPOVAL, 218 : & MELALV.MPOVAL, 219 : & MELALP.MPOVAL, 220 : & MELALY.MPOVAL 221 : -INC SMLENTI 222 : -INC SMELEME 223 : -INC SMLREEL 224 : POINTEUR MLRECV.MLREEL, MLRECP.MLREEL 225 : -INC SMLMOTS 226 : C 227 : C**** Segments des fractions massiques gauche et droit 228 : C 229 : SEGMENT FRAMAS 230 : REAL*8 FRAMG(NESP), FRAMD(NESP) 231 : ENDSEGMENT 232 : C 233 : SEGACT MLMOT1 234 : C 235 : C 236 : C**** Initialisation des parametres d'erreur déjà faite, i.e. 237 : C 238 : C LOGNEG = .FALSE. 239 : C LOGBOR = .FALSE. 240 : C MESERR = ' ' 241 : C MOTERR(1:40) = MESERR(1:40) 242 : C VALER = 0.0D0 243 : C VAL1 = 0.0D0 244 : C VAL2 = 0.0D0 245 : C 246 : C AB 247 : C 2010/09/14: LOGI1 computed can destroy the symmetry! 248 : C LOGI1 = .FALSE. 249 : LOGI1 = .FALSE. 250 : C 251 : C**** KRIPAD pour la correspondance global/local de centre 252 : C 252 : CALL KRIPAD(ICEN,MLENT1) 253 : C 254 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements 255 : C 256 : C Si i est le numero global d'un noeud de ICEN, 257 : C MLENT1.LECT(i) contient sa position, i.e. 258 : C 259 : C I = numero global du noeud centre 260 : C MLENT1.LECT(i) = numero local du noeud centre 261 : C 262 : C MLENT1 déjà activé, i.e. 263 : C 264 : C SEGACT MLENT1 265 : C 266 : C**** Activation de CHPOINTs 267 : C 268 : C densité + grad + limiteur 269 : C vitesse + grad + limiteur 270 : C pression + grad + limiteur 271 : C fract.mass + grad + limiteur 272 : C gamma 273 : C cosinus directeurs des normales aux surface 274 : C 274 : CALL LICHT(IROC , MPROC , TYPE, IGEOM) 274 : CALL LICHT(IGRROC, MPGRR , TYPE, IGEOM) 274 : CALL LICHT(IVITC, MPVITC , TYPE, IGEOM) 274 : CALL LICHT(IGRVC, MPGRV , TYPE, IGEOM) 274 : CALL LICHT(IPC , MPPC , TYPE, IGEOM) 274 : CALL LICHT(IGRPC, MPGRP , TYPE, IGEOM) 274 : CALL LICHT(IYC , MPYC , TYPE, IGEOM) 274 : CALL LICHT(IGRYC, MPGRY , TYPE, IGEOM) 274 : CALL LICHT(IGAMC, MPGAMC , TYPE, IGEOM) 274 : CALL LICHT(INORM, MPNORM , TYPE, IGEOM) 275 : C 276 : C**** Les MPOVALs 'Prediction' 277 : C 278 : IF(LOGTEM)THEN 279 : SEGINI, MPROP = MPROC 280 : SEGINI, MPPP = MPPC 281 : SEGINI, MPVITP = MPVITC 282 : SEGINI, MPYP = MPYC 283 : ELSE 284 : MPROP = MPROC 285 : MPPP = MPPC 286 : MPVITP = MPVITC 287 : MPYP = MPYC 288 : ENDIF 289 : C 290 : C**** Les Limiteurs 291 : C 291 : CALL LICHT(IALROC, MELALR , TYPE, IGEOM) 291 : CALL LICHT(IALVC, MELALV , TYPE, IGEOM) 291 : CALL LICHT(IALPC, MELALP , TYPE, IGEOM) 291 : CALL LICHT(IALYC, MELALY , TYPE, IGEOM) 292 : C 293 : C 294 : C**** Les MPOVAL sont déjà activés i.e.: 295 : C 296 : C SEGACT MPROC 297 : C SEGACT MPGRR 298 : C SEGACT MPIALR 299 : C SEGACT MPVITC 300 : C SEGACT MPGRV 301 : C SEGACT MPIALV 302 : C SEGACT MPPC 303 : C SEGACT MPGRP 304 : C SEGACT MPIALP 305 : C SEGACT MPYC 306 : C SEGACT MPGRY 307 : C SEGACT MPIALY 308 : C SEGACT MPGAMC 309 : C SEGACT MPNORM 310 : C 311 : C**** Le MELEME FACEL 312 : C 313 : IPT1 = IFACEL 314 : IPT2 = IFACE 315 : SEGACT IPT1 316 : SEGACT IPT2 317 : NFAC = IPT1.NUM(/2) 318 : C 319 : C**** Creation de MCHAMLs contenant les etats gauche et droite, 320 : C 321 : C i.e.: 322 : C 323 : C vitesse + cosinus directors du repere local 324 : C densité 325 : C pression 326 : C gamma 327 : C 328 : C**** Cosinus directors du repere local et vitesse 329 : C 330 : C Les cosinus directeurs 331 : C 332 : N1 = 2 333 : N3 = 6 334 : L1 = 28 335 : SEGINI MCHEL1 336 : IVITF = MCHEL1 337 : MCHEL1.TITCHE = 'U ' 338 : MCHEL1.IMACHE(1) = IFACE 339 : MCHEL1.IMACHE(2) = IFACEL 340 : MCHEL1.CONCHE(1) = ' (n,t) in (x,y) ' 341 : MCHEL1.CONCHE(2) = ' U in (n,t) ' 342 : C 343 : C**** Valeurs des cosinus definies par respect au repair global, i.e. 344 : C 345 : MCHEL1.INFCHE(1,1) = 2 346 : MCHEL1.INFCHE(1,3) = NIFOUR 347 : MCHEL1.INFCHE(1,4) = 0 348 : MCHEL1.INFCHE(1,5) = 0 349 : MCHEL1.INFCHE(1,6) = 0 350 : MCHEL1.IFOCHE = IFOUR 351 : C 352 : C**** Valeurs de vitesse definies par respect au repair local, i.e. 353 : C 354 : MCHEL1.INFCHE(2,1) = 1 355 : MCHEL1.INFCHE(2,3) = NIFOUR 356 : MCHEL1.INFCHE(2,4) = 0 357 : MCHEL1.INFCHE(2,5) = 0 358 : MCHEL1.INFCHE(2,6) = 0 359 : C 360 : C**** Le cosinus directeurs 361 : C 362 : N1PTEL = 1 363 : N1EL = NFAC 364 : N2PTEL = 0 365 : N2EL = 0 366 : C 367 : C**** MCHAML a N2 composantes: 368 : C 369 : C cosinus directeurs du repere local (n,t1) 370 : C 371 : C IDIM = 2 -> 4 composantes 372 : C 373 : N2 = 4 374 : SEGINI MCHAM1 375 : MCHEL1.ICHAML(1) = MCHAM1 376 : MCHAM1.NOMCHE(1) = 'NX ' 377 : MCHAM1.NOMCHE(2) = 'NY ' 378 : MCHAM1.NOMCHE(3) = 'TX ' 379 : MCHAM1.NOMCHE(4) = 'TY ' 380 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 381 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 382 : MCHAM1.TYPCHE(3) = 'REAL*8 ' 383 : MCHAM1.TYPCHE(4) = 'REAL*8 ' 384 : SEGINI MELVNX 385 : SEGINI MELVNY 386 : SEGINI MELT1X 387 : SEGINI MELT1Y 388 : MCHAM1.IELVAL(1) = MELVNX 389 : MCHAM1.IELVAL(2) = MELVNY 390 : MCHAM1.IELVAL(3) = MELT1X 391 : MCHAM1.IELVAL(4) = MELT1Y 392 : SEGDES MCHAM1 393 : C 394 : C**** Vitesse 395 : C 396 : N1EL = NFAC 397 : N1PTEL = 3 398 : N2EL = 0 399 : N2PTEL = 0 400 : C 401 : C**** MCHAML a N2 composantes: 402 : C 403 : C IDIM = 2 -> 2 composantes 404 : C 405 : N2 = 2 406 : SEGINI MCHAM1 407 : MCHEL1.ICHAML(2) = MCHAM1 408 : SEGDES MCHEL1 409 : MCHAM1.NOMCHE(1) = 'UN ' 410 : MCHAM1.NOMCHE(2) = 'UT ' 411 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 412 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 413 : SEGINI MELVUN 414 : SEGINI MELVUT 415 : MCHAM1.IELVAL(1) = MELVUN 416 : MCHAM1.IELVAL(2) = MELVUT 417 : SEGDES MCHAM1 418 : C 419 : C**** Densite 420 : C 421 : N1 = 1 422 : N3 = 6 423 : L1 = 15 424 : SEGINI MCHEL2 425 : IROF = MCHEL2 426 : MCHEL2.IMACHE(1) = IFACEL 427 : MCHEL2.TITCHE = 'RO ' 428 : MCHEL2.CONCHE(1) = ' ' 429 : C 430 : C**** Valeurs independente du repére, i.e. 431 : C 432 : MCHEL2.INFCHE(1,1) = 0 433 : MCHEL2.INFCHE(1,3) = NIFOUR 434 : MCHEL2.INFCHE(1,4) = 0 435 : MCHEL2.INFCHE(1,5) = 0 436 : MCHEL2.INFCHE(1,6) = 0 437 : MCHEL2.IFOCHE = IFOUR 438 : N2 = 1 439 : SEGINI MCHAM1 440 : MCHEL2.ICHAML(1) = MCHAM1 441 : SEGDES MCHEL2 442 : MCHAM1.NOMCHE(1) = 'SCAL ' 443 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 444 : SEGINI MELRO 445 : MCHAM1.IELVAL(1) = MELRO 446 : SEGDES MCHAM1 447 : C 448 : C**** Pression 449 : C 450 : MCHEL1 = IROF 451 : SEGINI, MCHEL2 = MCHEL1 452 : IPF = MCHEL2 453 : MCHEL2.TITCHE = 'P ' 454 : C 455 : C**** MCHAM1 = MCHAML de la densite 456 : C 457 : SEGINI, MCHAM2 = MCHAM1 458 : MCHEL2.ICHAML(1) = MCHAM2 459 : SEGDES MCHEL2 460 : SEGINI MELP 461 : MCHAM2.IELVAL(1) = MELP 462 : SEGDES MCHAM2 463 : C 464 : C**** Gamma 465 : C 466 : MCHEL1 = IROF 467 : SEGINI, MCHEL2 = MCHEL1 468 : IGAMF = MCHEL2 469 : MCHEL2.TITCHE = 'GAMMA ' 470 : C 471 : C**** MCHAM1 = MCHAML de la densite 472 : C 473 : SEGINI, MCHAM2 = MCHAM1 474 : MCHEL2.ICHAML(1) = MCHAM2 475 : SEGDES MCHEL2 476 : SEGINI MELGAM 477 : MCHAM2.IELVAL(1) = MELGAM 478 : SEGDES MCHAM2 479 : C 480 : C**** Le CHAMELEM des fractions massiques 481 : C 482 : NESP = MPYC.VPOCHA(/2) 483 : SEGINI FRAMAS 484 : MCHEL1 = IROF 485 : SEGINI, MCHEL2 = MCHEL1 486 : IYF = MCHEL2 487 : MCHEL2.TITCHE = 'Y ' 488 : N2 = NESP 489 : SEGINI MCHAMY 490 : MCHEL2.ICHAML(1) = MCHAMY 491 : SEGDES MCHEL2 492 : N1EL = NFAC 493 : N1PTEL = 3 494 : N2EL = 0 495 : N2PTEL = 0 496 : DO I1 = 1, NESP 497 : SEGINI MELVA1 498 : MCHAMY.IELVAL(I1) = MELVA1 499 : C CARCEL = MLMOT1.MOTS(1) 500 : C AB Error. It is corrected by the following line 501 : CARCEL = MLMOT1.MOTS(I1) 502 : TYPE = ' ' 503 : TYPE(1:4) = CARCEL(1:4) 504 : MCHAMY.NOMCHE(I1) = TYPE 505 : MCHAMY.TYPCHE(I1) = 'REAL*8 ' 506 : ENDDO 507 : C 508 : C 509 : C*********************************************************************** 510 : C********* PREDICTION ************************************************** 511 : C*********************************************************************** 512 : C 513 : C**** Prediction avec gradients limités 514 : C 515 : C 516 : IF(LOGTEM)THEN 517 : C 518 : IPT3 = ICEN 519 : SEGACT IPT3 520 : NCEN = IPT3.NUM(/2) 521 : DO NLCE = 1, NCEN 522 : ROG = MPROP.VPOCHA(NLCE,1) 523 : UXG = MPVITP.VPOCHA(NLCE,1) 524 : UYG = MPVITP.VPOCHA(NLCE,2) 525 : PG = MPPP.VPOCHA(NLCE,1) 526 : DROX = MPGRR.VPOCHA(NLCE,1)*MELALR.VPOCHA(NLCE,1) 527 : DROY = MPGRR.VPOCHA(NLCE,2)*MELALR.VPOCHA(NLCE,1) 528 : DUXX = MPGRV.VPOCHA(NLCE,1)*MELALV.VPOCHA(NLCE,1) 529 : DUXY = MPGRV.VPOCHA(NLCE,2)*MELALV.VPOCHA(NLCE,1) 530 : DUYX = MPGRV.VPOCHA(NLCE,3)*MELALV.VPOCHA(NLCE,2) 531 : DUYY = MPGRV.VPOCHA(NLCE,4)* MELALV.VPOCHA(NLCE,2) 532 : DPX = MPGRP.VPOCHA(NLCE,1)* MELALP.VPOCHA(NLCE,1) 533 : DPY = MPGRP.VPOCHA(NLCE,2)* MELALP.VPOCHA(NLCE,1) 534 : GAMG = MPGAMC.VPOCHA(NLCE,1) 535 : DRO = UXG * DROX + ROG * ( DUXX + DUYY ) 536 : & + UYG * DROY 537 : DUX = UXG * DUXX + DPX / ROG + UYG * DUXY 538 : DUY = UXG * DUYX + UYG * DUYY + DPY / ROG 539 : DP = GAMG * PG * (DUXX + DUYY) 540 : & + UXG * DPX + UYG * DPY 541 : C 542 : MPROP.VPOCHA(NLCE,1) = ROG - DELTAT * DRO 543 : MPVITP.VPOCHA(NLCE,1) = UXG - DELTAT * DUX 544 : MPVITP.VPOCHA(NLCE,2) = UYG - DELTAT * DUY 545 : MPPP.VPOCHA(NLCE,1) = PG - DELTAT * DP 546 : DO I1 = 1, NESP 547 : INDCEL = 2 * I1 - 1 548 : ALPHA = MELALY.VPOCHA(NLCE,I1) 549 : DYMAS = UXG * MPGRY.VPOCHA(NLCE,INDCEL) * ALPHA + 550 : & UYG * MPGRY.VPOCHA(NLCE,INDCEL+1) * ALPHA 551 : MPYP.VPOCHA(NLCE,I1) = MPYC.VPOCHA(NLCE,I1) - 552 : & DELTAT * DYMAS 553 : ENDDO 554 : ENDDO 555 : C 556 : ENDIF 557 : C 558 : C 559 : C*********************************************************************** 560 : C********* CORRECTION ************************************************** 561 : C*********************************************************************** 562 : C 563 : C**** Boucle sur le faces 564 : C 565 : IDIMP1 = IDIM + 1 566 : DO NLCF = 1, NFAC 567 : C 568 : C******* NLCF = numero local du centre de face 569 : C NGCF = numero global du centre de face 570 : C NGCEG = numero global du centre ELT "gauche" 571 : C NLCEG = numero local du centre ELT "gauche" 572 : C NGCED = numero global du centre ELT "droite" 573 : C NLCED = numero local du centre ELT "droite" 574 : C 575 : NGCEG = IPT1.NUM(1,NLCF) 576 : NGCF = IPT1.NUM(2,NLCF) 577 : NGCED = IPT1.NUM(3,NLCF) 578 : NLCEG = MLENT1.LECT(NGCEG) 579 : NLCED = MLENT1.LECT(NGCED) 580 : C 581 : C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF) 582 : C 583 : NGCF1 = IPT2.NUM(1,NLCF) 584 : IF( NGCF1 .NE. NGCF) THEN 585 : LOGAN = .TRUE. 586 : MESERR(1:40) = 'PRET, subroutine pre221.eso ' 587 : GOTO 9999 588 : ENDIF 589 : C 590 : C******* Cosinus directeurs des NORMALES aux faces 591 : C 592 : C On impose que les normales sont direct "Gauche" -> "Centre" 593 : C 594 : INDCEL = (NGCEG-1)*IDIMP1 595 : XG = XCOOR(INDCEL+1) 596 : YG = XCOOR(INDCEL+2) 597 : INDCEL = (NGCF-1)*IDIMP1 598 : XC = XCOOR(INDCEL + 1) 599 : YC = XCOOR(INDCEL + 2) 600 : INDCEL = (NGCED-1)*IDIMP1 601 : XD = XCOOR(INDCEL+1) 602 : YD = XCOOR(INDCEL+2) 603 : DXG = XC - XG 604 : DYG = YC - YG 605 : DXD = XC - XD 606 : DYD = YC - YD 607 : C 608 : C******* On calcule le sign du pruduit scalare 609 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 610 : C 611 : CNX = MPNORM.VPOCHA(NLCF,1) 612 : CNY = MPNORM.VPOCHA(NLCF,2) 613 : ORIENT = CNX * DXG + CNY * DYG 614 : ORIENT = SIGN(1.0D0,ORIENT) 615 : IF(ORIENT .NE. 1.0D0)THEN 616 : LOGAN = .TRUE. 617 : MESERR(1:30)= 618 : & 'PRET , subroutine pre221.eso. ' 619 : GOTO 9999 620 : ENDIF 621 : CNX = CNX * ORIENT 622 : CNY = CNY * ORIENT 623 : C 624 : C********** Cosinus directeurs de tangent 2D 625 : C 626 : CTX = -1.0D0 * CNY 627 : CTY = CNX 628 : C 629 : C 630 : C******* Les autres MELVALs 631 : C 632 : C 633 : C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0 634 : C GAMMA \in (1,3) 635 : C Si non il faut le faire, en utilisant LOGBOR, 636 : C LOGNEG, VALER, VAL1, VAL2 637 : C 638 : C 639 : C 640 : C******* NGCEG = NGCED -> Mur 641 : C 642 : IF( NGCEG .EQ. NGCED)THEN 643 : C 644 : C********** Sur le mur on fait de reconstruction sur l'etat gauche 645 : C 646 : C 647 : C********** Etat gauche 648 : C 649 : VALCEL = MPROP.VPOCHA(NLCEG, 1) 650 : ALCEL = MELALR.VPOCHA(NLCEG, 1) 651 : DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG + 652 : & MPGRR.VPOCHA(NLCEG, 2)*DYG 653 : ROG = VALCEL + ALCEL * DCEL 654 : C 655 : VALCEL = MPPP.VPOCHA(NLCEG, 1) 656 : ALCEL = MELALP.VPOCHA(NLCEG, 1) 657 : DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG + 658 : & MPGRP.VPOCHA(NLCEG, 2)*DYG 659 : PG = VALCEL + ALCEL * DCEL 660 : C 661 : LOGI2 = .FALSE. 662 : SUMY = 0.0D0 663 : CPCEL = 0.0D0 664 : CVCEL = 0.0D0 665 : DO I1 = 1, NESP 666 : INDCEL = 2 * I1 - 1 667 : VALCEL = MPYP.VPOCHA(NLCEG,I1) 668 : ALCEL = MELALY.VPOCHA(NLCEG, I1) 669 : DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG + 670 : & MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG 671 : ALCEL = VALCEL + ALCEL * DCEL 672 : SUMY = SUMY + ALCEL 673 : CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL 674 : CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL 675 : LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0) 676 : FRAMAS.FRAMG(I1) = ALCEL 677 : ENDDO 678 : LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0) 679 : CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1) 680 : CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1) 681 : GAMG = CPCEL / CVCEL 682 : C 683 : VALCEL = MPVITP.VPOCHA(NLCEG, 1) 684 : ALCEL = MELALV.VPOCHA(NLCEG, 1) 685 : DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG + 686 : & MPGRV.VPOCHA(NLCEG, 2)*DYG 687 : UXG = VALCEL + ALCEL * DCEL 688 : C 689 : VALCEL = MPVITP.VPOCHA(NLCEG, 2) 690 : ALCEL = MELALV.VPOCHA(NLCEG, 2) 691 : DCEL = MPGRV.VPOCHA(NLCEG, 3)*DXG + 692 : & MPGRV.VPOCHA(NLCEG, 4)*DYG 693 : UYG = VALCEL + ALCEL * DCEL 694 : C 695 : UNG = UXG * CNX + UYG * CNY 696 : UTG = UXG * CTX + UYG * CTY 697 : C 698 : C********** Si l'on fait pas de prediction, ce n'est pas necessaire de 699 : C controller la positivite' de la pression et de la densité; elle 700 : C est déjà garantie par la proprieté LED de limiteur; mais, vu 701 : C que le limiteur n'est pas calculé ici, mais dans un autre 702 : C operateur, on le fait 703 : C 704 : C LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2 705 : C 2010/09/14: AB this destroys the symmetry! 706 : C LOGI1 = .FALSE. 707 : C 708 : C 709 : IF(LOGI1)THEN 710 : C 711 : C************* Premier ordre en espace local 712 : C 713 : ROG = MPROC.VPOCHA(NLCEG,1) 714 : ROD = ROG 715 : PG = MPPC.VPOCHA(NLCEG,1) 716 : PD = PG 717 : UNG = MPVITC.VPOCHA(NLCEG,1)*CNX + 718 : & MPVITC.VPOCHA(NLCEG,2)*CNY 719 : UTG = MPVITC.VPOCHA(NLCEG,1)*CTX + 720 : & MPVITC.VPOCHA(NLCEG,2)*CTY 721 : UND = -1.0D0 * UNG 722 : UTD = UTG 723 : DO I1 = 1, NESP 724 : FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1) 725 : FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1) 726 : ENDDO 727 : GAMG = MPGAMC.VPOCHA(NLCEG,1) 728 : GAMD = GAMG 729 : ELSE 730 : C 731 : C********** Son etat droite 732 : C 733 : ROD = ROG 734 : PD = PG 735 : UND = -1.0D0 * UNG 736 : UTD = UTG 737 : DO I1 = 1, NESP 738 : FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1) 739 : ENDDO 740 : GAMD = GAMG 741 : ENDIF 742 : C 743 : C************* Fin cas mur 744 : C 745 : ELSE 746 : C 747 : C************* Etat gauche 748 : C 749 : VALCEL = MPROP.VPOCHA(NLCEG, 1) 750 : ALCEL = MELALR.VPOCHA(NLCEG, 1) 751 : DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG + 752 : & MPGRR.VPOCHA(NLCEG, 2)*DYG 753 : ROG = VALCEL + ALCEL * DCEL 754 : C 755 : VALCEL = MPPP.VPOCHA(NLCEG, 1) 756 : ALCEL = MELALP.VPOCHA(NLCEG, 1) 757 : DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG + 758 : & MPGRP.VPOCHA(NLCEG, 2)*DYG 759 : PG = VALCEL + ALCEL * DCEL 760 : C 761 : LOGI2 = .FALSE. 762 : SUMY = 0.0D0 763 : CPCEL = 0.0D0 764 : CVCEL = 0.0D0 765 : DO I1 = 1, NESP 766 : INDCEL = 2 * I1 - 1 767 : VALCEL = MPYP.VPOCHA(NLCEG,I1) 768 : ALCEL = MELALY.VPOCHA(NLCEG, I1) 769 : DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG + 770 : & MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG 771 : ALCEL = VALCEL + ALCEL * DCEL 772 : SUMY = SUMY + ALCEL 773 : CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL 774 : CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL 775 : LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0) 776 : FRAMAS.FRAMG(I1) = ALCEL 777 : ENDDO 778 : LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0) 779 : CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1) 780 : CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1) 781 : GAMG = CPCEL / CVCEL 782 : C 783 : VALCEL = MPVITP.VPOCHA(NLCEG, 1) 784 : ALCEL = MELALV.VPOCHA(NLCEG, 1) 785 : DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG + 786 : & MPGRV.VPOCHA(NLCEG, 2)*DYG 787 : UXG = VALCEL + ALCEL * DCEL 788 : C 789 : VALCEL = MPVITP.VPOCHA(NLCEG, 2) 790 : ALCEL = MELALV.VPOCHA(NLCEG, 2) 791 : DCEL = MPGRV.VPOCHA(NLCEG, 3)*DXG + 792 : & MPGRV.VPOCHA(NLCEG, 4)*DYG 793 : UYG = VALCEL + ALCEL * DCEL 794 : C 795 : UNG = UXG * CNX + UYG * CNY 796 : UTG = UXG * CTX + UYG * CTY 797 : C 798 : C********** Positivite 799 : C 800 : C LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2 801 : C 2010/09/14: AB this destroys the symmetry! 802 : C LOGI1 = .FALSE. 803 : C 804 : IF(LOGI1)THEN 805 : C 806 : C************* Premier ordre en espace local 807 : C 808 : ROG = MPROC.VPOCHA(NLCEG,1) 809 : PG = MPPC.VPOCHA(NLCEG,1) 810 : UNG = MPVITC.VPOCHA(NLCEG,1)*CNX + 811 : & MPVITC.VPOCHA(NLCEG,2)*CNY 812 : UTG = MPVITC.VPOCHA(NLCEG,1)*CTX + 813 : & MPVITC.VPOCHA(NLCEG,2)*CTY 814 : DO I1 = 1, NESP 815 : FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1) 816 : ENDDO 817 : GAMG = MPGAMC.VPOCHA(NLCEG,1) 818 : ENDIF 819 : C 820 : C********** Etat droite 821 : C 822 : VALCEL = MPROP.VPOCHA(NLCED, 1) 823 : ALCEL = MELALR.VPOCHA(NLCED, 1) 824 : DCEL = MPGRR.VPOCHA(NLCED, 1)*DXD + 825 : & MPGRR.VPOCHA(NLCED, 2)*DYD 826 : ROD = VALCEL + ALCEL * DCEL 827 : C 828 : VALCEL = MPPP.VPOCHA(NLCED, 1) 829 : ALCEL = MELALP.VPOCHA(NLCED, 1) 830 : DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD + 831 : & MPGRP.VPOCHA(NLCED, 2)*DYD 832 : PD = VALCEL + ALCEL * DCEL 833 : C 834 : LOGI2 = .FALSE. 835 : SUMY = 0.0D0 836 : CPCEL = 0.0D0 837 : CVCEL = 0.0D0 838 : DO I1 = 1, NESP 839 : INDCEL = 2 * I1 - 1 840 : VALCEL = MPYP.VPOCHA(NLCED,I1) 841 : ALCEL = MELALY.VPOCHA(NLCED, I1) 842 : DCEL = MPGRY.VPOCHA(NLCED, INDCEL)*DXD + 843 : & MPGRY.VPOCHA(NLCED,INDCEL + 1 )*DYD 844 : ALCEL = VALCEL + ALCEL * DCEL 845 : SUMY = SUMY + ALCEL 846 : CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL 847 : CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL 848 : LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0) 849 : FRAMAS.FRAMD(I1) = ALCEL 850 : ENDDO 851 : LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0) 852 : CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1) 853 : CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1) 854 : GAMD = CPCEL / CVCEL 855 : C 856 : VALCEL = MPVITP.VPOCHA(NLCED, 1) 857 : ALCEL = MELALV.VPOCHA(NLCED, 1) 858 : DCEL = MPGRV.VPOCHA(NLCED, 1)*DXD + 859 : & MPGRV.VPOCHA(NLCED, 2)*DYD 860 : UXD = VALCEL + ALCEL * DCEL 861 : C 862 : VALCEL = MPVITP.VPOCHA(NLCED, 2) 863 : ALCEL = MELALV.VPOCHA(NLCED, 2) 864 : DCEL = MPGRV.VPOCHA(NLCED, 3)*DXD + 865 : & MPGRV.VPOCHA(NLCED, 4)*DYD 866 : UYD = VALCEL + ALCEL * DCEL 867 : C 868 : UND = UXD * CNX + UYD * CNY 869 : UTD = UXD * CTX + UYD * CTY 870 : C 871 : C********** Positivite 872 : C 873 : C LOGI1 = (PD .LT. 0.0D0) .OR. (ROD .LT. 0.0D0) .OR. LOGI2 874 : C 2010/09/14: AB this destroys the symmetry! 875 : C LOGI1 = .FALSE. 876 : C 877 : IF(LOGI1)THEN 878 : C 879 : C************* Premier ordre en espace local 880 : C 881 : ROD = MPROC.VPOCHA(NLCED,1) 882 : PD = MPPC.VPOCHA(NLCED,1) 883 : UND = MPVITC.VPOCHA(NLCED,1)*CNX + 884 : & MPVITC.VPOCHA(NLCED,2)*CNY 885 : UTD = MPVITC.VPOCHA(NLCED,1)*CTX + 886 : & MPVITC.VPOCHA(NLCED,2)*CTY 887 : DO I1 = 1, NESP 888 : FRAMAS.FRAMD(I1) = MPYC.VPOCHA(NLCED,I1) 889 : ENDDO 890 : GAMD = MPGAMC.VPOCHA(NLCED,1) 891 : ENDIF 892 : ENDIF 893 : C 894 : C******** Les MELVALs 895 : C 896 : MELRO.VELCHE(1,NLCF) = ROG 897 : MELRO.VELCHE(3,NLCF) = ROD 898 : MELP.VELCHE(1,NLCF) = PG 899 : MELP.VELCHE(3,NLCF) = PD 900 : MELGAM.VELCHE(1,NLCF) = GAMG 901 : MELGAM.VELCHE(3,NLCF) = GAMD 902 : MELVUN.VELCHE(1,NLCF) = UNG 903 : MELVUN.VELCHE(3,NLCF) = UND 904 : MELVUT.VELCHE(1,NLCF) = UTG 905 : MELVUT.VELCHE(3,NLCF) = UTD 906 : MELVNX.VELCHE(1,NLCF) = CNX 907 : MELVNY.VELCHE(1,NLCF) = CNY 908 : MELT1X.VELCHE(1,NLCF) = CTX 909 : MELT1Y.VELCHE(1,NLCF) = CTY 910 : DO I1 = 1, NESP 911 : MELVA1 = MCHAMY.IELVAL(I1) 912 : MELVA1.VELCHE(1,NLCF) = FRAMAS.FRAMG(I1) 913 : MELVA1.VELCHE(3,NLCF) = FRAMAS.FRAMD(I1) 914 : ENDDO 915 : ENDDO 916 : C 917 : C**** Desactivation des SEGMENTs 918 : C 919 : SEGDES IPT1 920 : SEGDES IPT2 921 : C 922 : C**** Le MPOVALs 'Prediction' sont detruits (si existentes) 923 : C 924 : IF(LOGTEM)THEN 925 : SEGSUP MPROP 926 : SEGSUP MPVITP 927 : SEGSUP MPPP 928 : SEGSUP MPYP 929 : ENDIF 930 : C 931 : SEGDES MPROC 932 : SEGDES MPGRR 933 : SEGDES MELALR 934 : SEGDES MPVITC 935 : SEGDES MPGRV 936 : SEGDES MELALV 937 : SEGDES MPPC 938 : SEGDES MPGRP 939 : SEGDES MELALP 940 : SEGDES MPYC 941 : SEGDES MPGRY 942 : SEGDES MELALY 943 : SEGDES MPGAMC 944 : SEGDES MPNORM 945 : C 946 : SEGDES MELRO 947 : SEGDES MELP 948 : SEGDES MELGAM 949 : SEGDES MELVUN 950 : SEGDES MELVUT 951 : SEGDES MELVNX 952 : SEGDES MELVNY 953 : SEGDES MELT1X 954 : SEGDES MELT1Y 955 : C 956 : SEGDES MPYC 957 : DO I1 = 1, NESP 958 : MELVA1 = MCHAMY.IELVAL(I1) 959 : SEGDES MELVA1 960 : ENDDO 961 : SEGDES MCHAMY 962 : SEGSUP FRAMAS 963 : C 964 : SEGDES MLMOT1 965 : CC 966 : C**** Destruction du MELNTI correspondance local/global 967 : C 968 : SEGSUP MLENT1 969 : C 970 : 9999 CONTINUE 971 : C 972 : RETURN 973 : END 974 : 975 : 976 : 977 : 978 : 979 : 980 : 981 : 982 : 983 :
© Cast3M 2003 - All rights reserved.
Disclaimer