1 : C KFM12 SOURCE BECC 06/05/15 21:15:12 5442 2 : SUBROUTINE KFM12(IRES,IRN,IGN,IRETN,IGAMN,IVCO, 3 : & ICHPSU,ICHPDI,ICHPVO,INORM, 4 : & MELEMC,MELEMF,MELEFL,MELTFA,DTPS,DCFL, 5 : & MELLIM,ICHLIM,NJAC,ICOEF,ICOEV, 6 : & IDU,IPROBL) 7 : C************************************************************************ 8 : C 9 : C PROJET : CASTEM 2000 10 : C 11 : C NOM : KFM12 12 : C 13 : C DESCRIPTION : Voir aussi KFM1 14 : C 14 : C Cas deux dimensions, gaz "calorically perfect" 15 : C Méthode de relaxation de type SGS 16 : C 17 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 18 : C 19 : C AUTEUR : T. KLOCZKO DEN/DM2S/SFME/LTMF 20 : C 21 : C************************************************************************ 22 : C ENTREES 23 : C 24 : C 25 : C 1) Pointeurs des CHPOINTs 26 : C 27 : C IRES : CHPOINT 'CENTRE' contenant le residu 28 : C 29 : C IRN : CHPOINT 'CENTRE' contenant la masse volumique 30 : C 31 : C IGN : CHPOINT 'CENTRE' contenant la q.d.m. 32 : C 33 : C IRETN : CHPOINT 'CENTRE' contenant l'energie totale 34 : C 35 : C IGAMN : CHPOINT 'CENTRE' contenant le gamma 36 : C 37 : C IVCO : CHPOINT 'CENTRE' contenant le cut-off 38 : C 39 : C ICHLIM : CHPOINT conditions aux bords 40 : C 41 : C ICOEV : CHPOINT 'CENTRE' contenant le (gamma - 1) K / (R \rho) 42 : C coeff pour calculer le ray. spect. de la matrice 43 : C visc. 44 : C 45 : C 3) Pointeurs de CHPOINTs de la table DOMAINE 46 : C 47 : C ICHPSU : CHPOINT "FACE" contenant la surface des faces 48 : C 49 : C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum 50 : C de chaque element 51 : C 52 : C ICHPVO : CHPOINT "CENTRE" contenant le volume 53 : C de chaque element 54 : C 55 : C INORM : CHPOINT "CENTRE" contenant le normales aux faces 56 : C 57 : C 58 : C 4) Pointeurs de MELEME de la table DOMAINE 59 : C 60 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES 61 : C 62 : C MELEMF : MELEME 'FACE' du SPG des FACES 63 : C 64 : C MELEFL : MELEME 'FACEL' du connectivité FACES -> ELEM 65 : C 66 : C MELLIM : MELEME SPG conditions aux bords 67 : C 68 : C 5) 69 : C 70 : C DCFL = le double de la CFL 71 : C 72 : C DTPS = pas de temps physique 73 : C 74 : C NJAC = nombre d'iterations dans le jacobien 75 : C 76 : C 77 : C SORTIES 78 : C 79 : C IDU : resultat 80 : C 81 : C IPROBL : si > 0, il y a un probleme 82 : C 83 : C************************************************************************ 84 : C 85 : C HISTORIQUE (Anomalies et modifications éventuelles) 86 : C 87 : C HISTORIQUE : 2004 : creation 88 : C 89 : C 90 : C************************************************************************ 91 : C 92 : C 93 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0 94 : C GAMMA \in (1,3) 95 : C Y \in (0,1) 96 : C Si non il faut le faire!!! 97 : C 98 : C************************************************************************ 99 : C 100 : C 101 : C**** Variables de COOPTIO 102 : C 103 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 104 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 105 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 106 : C & ,IECHO, IIMPI, IOSPI 107 : C & ,IDIM 108 : C & ,MCOORD 109 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 110 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 111 : C & ,NORINC,NORVAL,NORIND,NORVAD 112 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI, IREFOR, ISAFOR 113 : C 114 : IMPLICIT INTEGER(I-N) 115 : INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPDI,ICHPVO,INORM 116 : & ,MELEMF,IDU,NFAC, NFACE 117 : & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCEG, NLCED 118 : & ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP 119 : & ,I1,I2, IBLJAC, IPROBL, IVCO, MELTFA 120 : & ,NSOUPO,IFACT,NELT, ISOUP, IELT, IFAC, IFACG, JG 121 : & ,ISOUP1, ISOUP2, ICOEF,ICOEV 122 : & , ID, IELT1, IELT2, IFAC1, IFAC2 123 : REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, RECG, RECD, PG, VOLG 124 : & , DIAMG, ROD, RUXD, RUYD, UND, RETD, GAMD, PD 125 : & , CNX, CNY, DCFL, DTPS 126 : & , SURF, SIGMAF 127 : & , UNSDT, UNSDTF, UXD, UYD 128 : & , FR, FRUX, FRUY, FRET, COEF 129 : & , CTX, CTY, RHTD, UTD, UTG, UXG, UYG, QG, CG, UCOG 130 : & , ROF, PF, UNF, GAMF, UCOF, CF, PHI2, COEF1 131 : & , QN 132 : C 133 : CHARACTER*8 TYPE 134 : LOGICAL LOGINV, LOGJAC 135 : C 136 : C**** LES INCLUDES 137 : C 138 : -INC CCOPTIO 139 : -INC SMCHPOI 140 : POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL 141 : & , MPODU.MPOVAL, MPOVOL.MPOVAL 142 : & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL 143 : & , MPNORM.MPOVAL, MPLIM.MPOVAL, MPORES.MPOVAL, MPOCO.MPOVAL 144 : & , MPOCO1.MPOVAL, MELEFL.MELEME, MELEMC.MELEME 145 : & , MPCOEV.MPOVAL 146 : -INC SMCOORD 147 : -INC SMELEME 148 : -INC SMLMOTS 149 : -INC SMLENTI 150 : POINTEUR MLELIM.MLENTI, MLESUP.MLENTI 151 : C 152 : C**** POINTEUR vecteurs 153 : C 154 : SEGMENT VEC 155 : REAL*8 VALVEC(I1) 156 : ENDSEGMENT 157 : POINTEUR SIGMA0.VEC, AN0.VEC, DN0.VEC, D1.VEC 158 : & , USDTI.VEC, PHI2V.VEC, SIGMAV.VEC 159 : C 160 : C**** POINTEUR matriciels 161 : C 162 : SEGMENT MAT 163 : REAL*8 VAL(I1,I2) 164 : ENDSEGMENT 165 : POINTEUR CONS0.MAT, CONS1.MAT, BN0.MAT, BN1.MAT 166 : & , CN0.MAT, CN1.MAT, DU.MAT 167 : & , RES.MAT 168 : & , Q1.MAT, Q2.MAT 169 : C 170 : C 171 : C**** IPROBL: si different de 0, un probleme a été detecté 172 : C 173 : IPROBL=0 174 : C 175 : C**** Initialisation des MLENTI des conditions aux limites 176 : C 177 : C 177 : CALL KRIPAD(MELLIM,MLELIM) 178 : C SEGINI MLELIM 179 : C 180 : C**** Initialisation des MELEMEs 181 : C 182 : C 'CENTRE', 'FACEL' 183 : C 184 : SEGACT MELEFL 185 : NFAC = MELEFL.NUM(/2) 186 : C 187 : C**** KRIPAD pour la correspondance global/local de centre 188 : C 188 : CALL KRIPAD(MELEMC,MLENT1) 189 : C 190 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements 191 : C 192 : C Si i est le numero global d'un noeud de ICEN, 193 : C MLENT1.LECT(i) contient sa position, i.e. 194 : C 195 : C I = numero global du noeud centre 196 : C MLENT1.LECT(i) = numero local du noeud centre 197 : C 198 : C MLENT1 déjà activé, i.e. 199 : C 200 : C SEGINI MLENT1 201 : C 202 : C 203 : C**** KRIPAD pour la correspondance global/local de 'FACE' 204 : C 204 : CALL KRIPAD(MELEMF,MLENT2) 205 : C SEGINI MLENT2 206 : C 207 : C 208 : C**** CHPOINTs de la table DOMAINE 209 : C 209 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF) 209 : CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC) 209 : CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC) 209 : CALL LICHT(INORM,MPNORM,TYPE,IGEOMC) 210 : C 211 : C**** LICHT active les MPOVALs en *MOD 212 : C 213 : C i.e. 214 : C 215 : C SEGACT MPOVSU*MOD 216 : C SEGACT MPOVOL*MOD 217 : C SEGACT MPOVDI*MOD 218 : C SEGACT MPNORM*MOD 219 : C 219 : CALL LICHT(IDU,MPODU,TYPE,IGEOMC) 220 : C SEGACT MPODU*MOD 221 : C 222 : C USDTI initialisé a zero; USDTI = 1 / DT 223 : C 224 : NCEN=MPOVOL.VPOCHA(/1) 225 : I1=NCEN 226 : SEGINI USDTI 227 : C 227 : CALL LICHT(IRES,MPORES,TYPE,IGEOMC) 227 : CALL LICHT(IRN,MPRN,TYPE,IGEOMC) 227 : CALL LICHT(IGN,MPGN,TYPE,IGEOMC) 227 : CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC) 227 : CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC) 227 : CALL LICHT(IVCO,MPOCO,TYPE,IGEOMC) 227 : CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM) 227 : CALL LICHT(ICOEV,MPCOEV,TYPE,IGEOMC) 228 : C 229 : C SEGACT MPORES*MOD 230 : C SEGACT MPRN*MOD 231 : C SEGACT MPGN*MOD 232 : C SEGACT MPRETN*MOD 233 : C SEGACT MPGAMN*MOD 234 : C SEGACT MPOCO 235 : C SEGACT MPLIM*MOD 236 : C SEGACT MPCOEV*MOD 237 : C 238 : SEGINI, MPOCO1=MPOCO 239 : C 240 : C**** MPOCO1 is the "corrected cut-off" speed. 241 : C It is bigger than the given cut-off and the local 242 : C speed. 243 : C 244 : IF(IGEOMF .NE. MELEMF)THEN 245 : WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso' 246 : WRITE(IOIMP,*) 'Probleme de SPG' 247 : C 21 2 248 : C Données incompatibles 248 : CALL ERREUR(21) 249 : GOTO 9999 250 : ENDIF 251 : IF(IGEOMC .NE. MELEMC)THEN 252 : WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso' 253 : WRITE(IOIMP,*) 'Probleme de SPG' 254 : C 21 2 255 : C Données incompatibles 255 : CALL ERREUR(21) 256 : GOTO 9999 257 : ENDIF 258 : C 259 : C**** On initialise GAMMA=I+(1-phi2)/phi2*Q1.(Q2)^t 260 : C 261 : I1=4 262 : I2=NCEN 263 : SEGINI Q1 264 : SEGINI Q2 265 : C 266 : C**** On initialise le segment de travail D1 267 : C 268 : I1=4 269 : SEGINI D1 270 : C 271 : C**** On initialise AN0, DN0, PHI2V 272 : C 273 : I1=NCEN 274 : SEGINI AN0 275 : SEGINI DN0 276 : SEGINI PHI2V 277 : C 278 : C**** On initialise BN0 279 : C 280 : I1=4 281 : I2=NCEN 282 : SEGINI BN0 283 : C 284 : C**** On initialise BN1 285 : C 286 : I1=4 287 : I2=NCEN 288 : SEGINI BN1 289 : C 290 : C**** On initialise CN0 291 : C 292 : I1=4 293 : I2=NCEN 294 : SEGINI CN0 295 : C 296 : C**** On initialise CN1 297 : C 298 : I1=4 299 : I2=NCEN 300 : SEGINI CN1 301 : C 302 : C**** On initialise SIGMA0, SIGMAV 303 : C 304 : I1=NFAC 305 : SEGINI SIGMA0 306 : SEGINI SIGMAV 307 : C 308 : C**** RES 309 : C 310 : I1=4 311 : I2=NCEN 312 : SEGINI RES 313 : DO ICEN=1,NCEN,1 314 : DO ICOMP=1,4,1 315 : RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP) 316 : ENDDO 317 : ENDDO 318 : SEGDES MPORES 319 : C 320 : C***** DU (increment des variables conservatives) 321 : C 322 : I1=4 323 : I2=NCEN 324 : SEGINI DU 325 : C 326 : C**** Les variables conservatives 327 : C 328 : C CONS1.VAL contient les variables conservatives + gamma 329 : C i.e. rho(i) = CONS.VAL(1,i) 330 : C rhoux(i) = CONS.VAL(2,i) 331 : C rhouy(i) = CONS.VAL(3,i) 332 : C rhoet(i) = CONS.VAL(4,i) 333 : C gam(i) = CONS.VAL(5,i) 334 : C 335 : I1=5 336 : I2=NCEN 337 : SEGINI CONS0 338 : DO ICEN=1,NCEN,1 339 : CONS0.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1) 340 : CONS0.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1) 341 : CONS0.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2) 342 : CONS0.VAL(4,ICEN)=MPRETN.VPOCHA(ICEN,1) 343 : CONS0.VAL(5,ICEN)=MPGAMN.VPOCHA(ICEN,1) 344 : ENDDO 345 : SEGDES MPRN 346 : SEGDES MPRETN 347 : SEGDES MPGN 348 : SEGDES MPGAMN 349 : C 350 : SEGINI, CONS1=CONS0 351 : C 352 : C**** Storage of l-1 values of F at the l Jacobi iteration 353 : C 354 : IPT1=MELTFA 355 : SEGACT IPT1 356 : NSOUPO= IPT1.LISOUS(/1) 357 : NSOUPO=MAX(NSOUPO,1) 358 : JG=NSOUPO 359 : SEGINI MLESUP 360 : IFACT=0 361 : IF(NSOUPO .EQ. 1)THEN 362 : MLESUP.LECT(1)=IPT1 363 : NFACE=IPT1.NUM(/1) 364 : NELT=IPT1.NUM(/2) 365 : IFACT=IFACT+(NFACE*NELT) 366 : ELSE 367 : DO ISOUP = 1, NSOUPO, 1 368 : IPT2=IPT1.LISOUS(ISOUP) 369 : MLESUP.LECT(ISOUP)=IPT2 370 : SEGACT IPT2 371 : NFACE=IPT2.NUM(/1) 372 : NELT=IPT2.NUM(/2) 373 : IFACT=IFACT+(NFACE*NELT) 374 : ENDDO 375 : ENDIF 376 : C 377 : C*********************************************************************** 378 : C**** FIRST JACOBI ITERATION ******************************************* 379 : C*********************************************************************** 380 : C 381 : C 382 : C**** LOOP on ELTFA to compute: 383 : C 384 : C AN0_i = 1/(2 \Omega_i) \sum_j \sigma_i,j S_j 385 : C SIGMA0_j 386 : C The cut-off in each cell 387 : C The reciprocal of the dual time step USDTI 388 : C 389 : ICEN=0 390 : IFACG=0 391 : DO ISOUP=1,NSOUPO,1 392 : IPT2=MLESUP.LECT(ISOUP) 393 : NFACE=IPT2.NUM(/1) 394 : NELT=IPT2.NUM(/2) 395 : DO IELT=1,NELT,1 396 : ICEN=ICEN+1 397 : DO IFAC=1,NFACE,1 398 : C IFACG increase with IFAC,IELT,ISOUP 399 : IFACG = IFACG+1 400 : NGCF = IPT2.NUM(IFAC,IELT) 401 : NLCF = MLENT2.LECT(NGCF) 402 : C 403 : C************* NLCF = numero local du centre de facel 404 : C NGCF = numero global du centre de facel 405 : C NGCEG = numero global du centre ELT "gauche" 406 : C NLCEG = numero local du centre ELT "gauche" 407 : C NGCED = numero global du centre ELT "droite" 408 : C NLCED = numero local du centre ELT "droite" 409 : C 410 : NGCEG = MELEFL.NUM(1,NLCF) 411 : NGCED = MELEFL.NUM(3,NLCF) 412 : NLCEG = MLENT1.LECT(NGCEG) 413 : NLCED = MLENT1.LECT(NGCED) 414 : C 415 : CNX = MPNORM.VPOCHA(NLCF,1) 416 : CNY = MPNORM.VPOCHA(NLCF,2) 417 : CTX = -1.0D0 * CNY 418 : CTY = CNX 419 : C 420 : C************* Left or right? 421 : C 422 : IF(NLCEG .EQ. ICEN)THEN 423 : LOGINV=.FALSE. 424 : ELSEIF(NLCED .EQ. ICEN)THEN 425 : LOGINV=.TRUE. 426 : ELSE 427 : WRITE(IOIMP,*) 'Anomalie dans kfm13.eso' 428 : WRITE(IOIMP,*) 'Probleme de SPG' 429 : C 21 2 430 : C Données incompatibles 430 : CALL ERREUR(21) 431 : GOTO 9999 432 : ENDIF 433 : C************* If right, inversion 434 : IF(LOGINV) THEN 435 : NLCED=NLCEG 436 : NLCEG=ICEN 437 : CNX = -1.0D0*CNX 438 : CNY = -1.0D0*CNY 439 : CTX = -1.0D0*CTX 440 : CTY = -1.0D0*CTY 441 : ENDIF 442 : C 443 : SURF = MPOVSU.VPOCHA(NLCF,1) 444 : C 445 : C************* Récupération des Etats "gauche" et "droite" 446 : C 447 : C 448 : ROG = CONS0.VAL(1,NLCEG) 449 : RUXG = CONS0.VAL(2,NLCEG) 450 : RUYG = CONS0.VAL(3,NLCEG) 451 : UXG = RUXG / ROG 452 : UYG = RUYG / ROG 453 : UNG = (UXG * CNX) + (UYG * CNY) 454 : UTG = (UXG * CTX) + (UYG * CTY) 455 : RETG = CONS0.VAL(4,NLCEG) 456 : GAMG = CONS0.VAL(5,NLCEG) 457 : RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG)) 458 : RECG = RECG / ROG 459 : PG = (GAMG - 1.0D0) * (RETG - RECG) 460 : VOLG = MPOVOL.VPOCHA(NLCEG,1) 461 : DIAMG = MPOVDI.VPOCHA(NLCEG,1) 462 : C 463 : ROD = CONS0.VAL(1,NLCED) 464 : RUXD = CONS0.VAL(2,NLCED) 465 : RUYD = CONS0.VAL(3,NLCED) 466 : UXD = RUXD / ROD 467 : UYD = RUYD / ROD 468 : RETD = CONS0.VAL(4,NLCED) 469 : GAMD = CONS0.VAL(5,NLCED) 470 : RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD)) 471 : RECD = RECD / ROD 472 : PD = (GAMD - 1.0D0) * (RETD - RECD) 473 : UND = (UXD * CNX) + (UYD * CNY) 474 : UTD = (UXD * CTX) + (UYD * CTY) 475 : C 476 : IF(NLCEG .EQ. NLCED)THEN 477 : C Murs ou condition limite 478 : NLLIM=MLELIM.LECT(NGCF) 479 : IF(NLLIM .EQ.0)THEN 480 : C Mur 481 : UND = -1.0D0 * UNG 482 : UTD = UTG 483 : UXD = UND * CNX + UTD * CTX 484 : UYD = UND * CNY + UTD * CTY 485 : RUXD = ROD*UXD 486 : RUYD = ROD*UYD 487 : ELSE 488 : ROD = MPLIM.VPOCHA(NLLIM,1) 489 : UXD = MPLIM.VPOCHA(NLLIM,2) 490 : UYD = MPLIM.VPOCHA(NLLIM,3) 491 : RUXD = ROD*UXD 492 : RUYD = ROD*UYD 493 : PD = MPLIM.VPOCHA(NLLIM,4) 494 : GAMD = GAMG 495 : UND = (UXD * CNX) + (UYD * CNY) 496 : UTD = (UXD * CTX) + (UYD * CTY) 497 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD) 498 : & +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD))) 499 : ENDIF 500 : ENDIF 501 : C 502 : C************* We check that density and pressure are positive 503 : C 504 : IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR. 505 : & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN 506 : WRITE(IOIMP,*) 'FMM' 507 : WRITE(IOIMP,*) 'TEST' 508 : C 22 0 509 : C**************** Opération malvenue. Résultat douteux 510 : C 510 : CALL ERREUR(22) 511 : IPROBL = 1 512 : GOTO 9998 513 : ENDIF 514 : C 515 : C************* We compute the interfacial state 516 : C 517 : ROF = 0.5D0*(ROG+ROD) 518 : PF = 0.5D0*(PG+PD) 519 : UNF = 0.5D0*(UNG+UND) 520 : GAMF = 0.5D0*(GAMG+GAMD) 521 : C 522 : C************* We compute SIGMA and the corrected cut-off speed 523 : C 524 : UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1)) 525 : UCOF=MAX(UCOF,(((UNG*UNG)+(UTG*UTG))**0.5D0)) 526 : UCOF=MAX(UCOF,(((UND*UND)+(UTD*UTD))**0.5D0)) 527 : IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1)) 528 : & MPOCO1.VPOCHA(NLCEG,1)= UCOF 529 : C 530 : QN=ABS(UNF) 531 : CF=(GAMF*PF/ROF)**0.5D0 532 : IF(UCOF .GE. CF)THEN 533 : PHI2 = 1.0D0 534 : ELSEIF(QN .GT. CF)THEN 535 : PHI2 = 1.0D0 536 : ELSEIF(QN .LT. UCOF)THEN 537 : PHI2 = UCOF / CF 538 : ELSE 539 : PHI2 = QN / CF 540 : ENDIF 541 : PHI2=PHI2*PHI2 542 : C 543 : SIGMAF=(1.0D0-PHI2)*QN 544 : SIGMAF=SIGMAF*SIGMAF 545 : SIGMAF=SIGMAF+(4.0D0*PHI2*CF*CF) 546 : SIGMAF=SIGMAF**0.5D0 547 : SIGMAF=SIGMAF+((1.0D0+PHI2)*QN) 548 : SIGMAF=0.5D0*SIGMAF 549 : C 550 : SIGMA0.VALVEC(NLCF)=SIGMAF 551 : C 552 : C************* We compute AN0 553 : C 554 : AN0.VALVEC(NLCEG)=AN0.VALVEC(NLCEG) 555 : & +(0.5D0*SIGMAF*SURF/VOLG) 556 : C 557 : C************* We compute the dual time step 558 : C 559 : UNSDT=USDTI.VALVEC(NLCEG) 560 : UNSDTF = SIGMAF / DIAMG 561 : IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG)=UNSDTF 562 : C 563 : C************* We compute SIGMAV 564 : C 565 : C COEF=FG 566 : COEF=(MCOORD.XCOOR(((NGCEG-1)*3)+1)- 567 : & MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX 568 : COEF=COEF+ 569 : & ((MCOORD.XCOOR(((NGCEG-1)*3)+2)- 570 : & MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY) 571 : C COEF1=FD 572 : COEF1=(MCOORD.XCOOR(((NGCED-1)*3)+1)- 573 : & MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX 574 : COEF1=COEF1+ 575 : & ((MCOORD.XCOOR(((NGCED-1)*3)+2)- 576 : & MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY) 577 : C COEF=|FG|+|FD| 578 : COEF=ABS(COEF)+ABS(COEF1) 579 : SIGMAF=0.5D0*(MPCOEV.VPOCHA(NLCEG,1)+ 580 : & MPCOEV.VPOCHA(NLCED,1))/COEF 581 : SIGMAV.VALVEC(NLCF)=SIGMAF 582 : C 583 : C************* We compute the contribution of SIGMAV to 584 : C the diagonal part of the matrix to invert 585 : C 586 : C In the given report, DNO is a_i^0, i.e. 587 : C DN0 = AN0 (previously defined) + 588 : C + (\frac{1}{{\Omega_i}} 589 : C \sum_j \sigma_{V,i,j}^{0} S_j) 590 : C + 1 / \Delta t^n 591 : C 592 : COEF=SURF/VOLG 593 : DN0.VALVEC(NLCEG)=DN0.VALVEC(NLCEG)+COEF*SIGMAF 594 : C 595 : C 596 : C************* We compute the flux 597 : C FLU.VALVEC(1:4,IFAC) flux in IFAC face 598 : C 1 -> mass 599 : C 2 -> momentum (x) 600 : C 3 -> momentum (y) 601 : C 4 -> energy 602 : C 603 : RHTD = RETD+PD 604 : FR = (ROD*UND) 605 : FRUX = (ROD*UND*UXD)+(PD*CNX) 606 : FRUY = (ROD*UND*UYD)+(PD*CNY) 607 : FRET = (UND*RHTD) 608 : C 609 : COEF = SURF/(2.0D0*VOLG) 610 : C 611 : BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG) 612 : & + FR * COEF 613 : BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG) 614 : & + FRUX * COEF 615 : BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG) 616 : & + FRUY * COEF 617 : BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG) 618 : & + FRET * COEF 619 : C 620 : C************* Viscous contribution 621 : C 622 : COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG 623 : C 624 : BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG) 625 : & - (COEF * ROD) 626 : BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG) 627 : & - (COEF * RUXD) 628 : BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG) 629 : & - (COEF * RUYD) 630 : BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG) 631 : & - (COEF * RETD) 632 : C 633 : C************* On calcule CN0 634 : C 635 : C CN0 in the referred report is given by 636 : C \sum_j 637 : C \left\{ 638 : C \frac{S_j}{2 \Omega_i} ~ 639 : C \left( 640 : C \sigma_{i,j}^{0} 641 : C \mathbf{U}_{j,R}^{n+1,m} 642 : C \right) 643 : C \right\} 644 : C 645 : C 646 : COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG) 647 : C 648 : CN0.VAL(1,NLCEG) = CN0.VAL(1,NLCEG) 649 : & + (COEF * ROD) 650 : CN0.VAL(2,NLCEG) = CN0.VAL(2,NLCEG) 651 : & + (COEF * RUXD) 652 : CN0.VAL(3,NLCEG) = CN0.VAL(3,NLCEG) 653 : & + (COEF * RUYD) 654 : CN0.VAL(4,NLCEG) = CN0.VAL(4,NLCEG) 655 : & + (COEF * RETD) 656 : 657 : ENDDO 658 : C 659 : C********** End of loop on the faces of the element 660 : C 661 : C N.B.: ICEN = NLCEG 662 : C 663 : C 664 : C Computation of the low-Mach preconditioning Matrix 665 : C 666 : CG=(GAMG*PG/ROG)**0.5D0 667 : UCOG=MPOCO1.VPOCHA(ICEN,1) 668 : QG=2.0D0*RECG/ROG 669 : QG=QG**0.5D0 670 : C 671 : IF(UCOG .GE. CG)THEN 672 : PHI2 = 1.0D0 673 : ELSEIF(QG .GT. CG)THEN 674 : PHI2 = 1.0D0 675 : ELSEIF(QG .LT. UCOG)THEN 676 : PHI2 = UCOG / CG 677 : ELSE 678 : PHI2 = QG / CG 679 : ENDIF 680 : PHI2=PHI2*PHI2 681 : PHI2V.VALVEC(ICEN)=PHI2 682 : C 683 : Q2.VAL(1,ICEN)=(0.5D0*QG*QG) 684 : Q2.VAL(2,ICEN)=-1.0D0*RUXG/ROG 685 : Q2.VAL(3,ICEN)=-1.0D0*RUYG/ROG 686 : Q2.VAL(4,ICEN)=1.0D0 687 : C 688 : COEF1=(GAMG-1.0D0)/(CG*CG) 689 : COEF=(1.0D0/COEF1)+(0.5D0*QG*QG) 690 : C 691 : Q1.VAL(1,ICEN)=COEF1 692 : Q1.VAL(2,ICEN)=COEF1*RUXG/ROG 693 : Q1.VAL(3,ICEN)=COEF1*RUYG/ROG 694 : Q1.VAL(4,ICEN)=COEF1*COEF 695 : C 696 : C 697 : C Computation of the diagonal coefficients 698 : C 699 : AN0.VALVEC(ICEN)=AN0.VALVEC(ICEN) 700 : C & +(USDTI.VALVEC(ICEN)*2.0D0/DCFL) 701 : C 702 : DN0.VALVEC(ICEN)=DN0.VALVEC(ICEN)+(1.0D0/DTPS) 703 : & +(USDTI.VALVEC(ICEN)*2.0D0/DCFL) 704 : C 705 : ENDDO 706 : ENDDO 707 : C 708 : C*********************************************************************** 709 : C**** END FIRST JACOBI ITERATION *************************************** 710 : C*********************************************************************** 711 : C 712 : C*********************************************************************** 713 : C**** JACOBI LOOP ****************************************************** 714 : C*********************************************************************** 715 : C 716 : C 717 : LOGJAC=.TRUE. 718 : C 719 : DO IBLJAC=1,NJAC,1 720 : C 721 : C 722 : C Si ICOEF=2 (i.e. si on a "LJACOF" (voir kfm1.eso)) 723 : C alors relaxation SGS avec uniquement des balayages "forward" 724 : C 725 : C Si ICOEF=3 (i.e. si on a "LJACOB") 726 : C alors relaxation SGS avec uniquement des balayages "backward" 727 : C 728 : C Sinon (i.e. si on a "LJACOFB") 729 : C alors relaxation SGS avec des balayages "forward" et "backward" 730 : C 731 : IF(ICOEF .EQ. 2)THEN 732 : LOGJAC=.TRUE. 733 : ELSEIF(ICOEF .EQ. 3)THEN 734 : LOGJAC=.FALSE. 735 : ELSE 736 : LOGINV=LOGJAC 737 : LOGJAC=(IBLJAC / 2 * 2) .EQ. IBLJAC 738 : IF(LOGJAC)THEN 739 : LOGJAC= .NOT. LOGINV 740 : ELSE 741 : LOGJAC= LOGINV 742 : ENDIF 743 : ENDIF 744 : C 745 : IF(LOGJAC)THEN 746 : ICEN=0 747 : IFACG=0 748 : ID=1 749 : ISOUP1=1 750 : ISOUP2=NSOUPO 751 : ELSE 752 : ICEN=NCEN+1 753 : IFACG=IFACT+1 754 : ID=-1 755 : ISOUP1=NSOUPO 756 : ISOUP2=1 757 : ENDIF 758 : C 759 : DO ISOUP=ISOUP1,ISOUP2,ID 760 : IPT2=MLESUP.LECT(ISOUP) 761 : NFACE=IPT2.NUM(/1) 762 : NELT=IPT2.NUM(/2) 763 : C 764 : IF(LOGJAC)THEN 765 : IELT1=1 766 : IELT2=NELT 767 : IFAC1=1 768 : IFAC2=NFACE 769 : ELSE 770 : IELT2=1 771 : IELT1=NELT 772 : IFAC2=1 773 : IFAC1=NFACE 774 : ENDIF 775 : C 776 : DO IELT=IELT1,IELT2,ID 777 : ICEN=ICEN+ID 778 : DO IFAC=IFAC1,IFAC2,ID 779 : C IFACG increase with IFAC,IELT,ISOUP 780 : IFACG = IFACG+ID 781 : NGCF = IPT2.NUM(IFAC,IELT) 782 : NLCF = MLENT2.LECT(NGCF) 783 : NGCEG = MELEFL.NUM(1,NLCF) 784 : NGCED = MELEFL.NUM(3,NLCF) 785 : NLCEG = MLENT1.LECT(NGCEG) 786 : NLCED = MLENT1.LECT(NGCED) 787 : C 788 : CNX = MPNORM.VPOCHA(NLCF,1) 789 : CNY = MPNORM.VPOCHA(NLCF,2) 790 : CTX = -1.0D0 * CNY 791 : CTY = CNX 792 : C 793 : C**************** Left or right? 794 : C 795 : IF(NLCEG .EQ. ICEN)THEN 796 : LOGINV=.FALSE. 797 : ELSEIF(NLCED .EQ. ICEN)THEN 798 : LOGINV=.TRUE. 799 : ELSE 800 : WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso' 801 : WRITE(IOIMP,*) 'Probleme de SPG' 802 : C 21 2 803 : C Données incompatibles 803 : CALL ERREUR(21) 804 : GOTO 9999 805 : ENDIF 806 : C**************** If right, inversion 807 : IF(LOGINV) THEN 808 : NLCED=NLCEG 809 : NLCEG=ICEN 810 : CNX = -1.0D0*CNX 811 : CNY = -1.0D0*CNY 812 : CTX = -1.0D0*CTX 813 : CTY = -1.0D0*CTY 814 : ENDIF 815 : C 816 : SURF = MPOVSU.VPOCHA(NLCF,1) 817 : SIGMAF = SIGMA0.VALVEC(NLCF) 818 : C 819 : C**************** Recuperation des Etats "gauche" et "droite" 820 : C 821 : ROG = CONS1.VAL(1,NLCEG) 822 : RUXG = CONS1.VAL(2,NLCEG) 823 : RUYG = CONS1.VAL(3,NLCEG) 824 : ROG = CONS1.VAL(1,NLCEG) 825 : RUXG = CONS1.VAL(2,NLCEG) 826 : RUYG = CONS1.VAL(3,NLCEG) 827 : UXG = RUXG / ROG 828 : UYG = RUYG / ROG 829 : UNG = (UXG * CNX) + (UYG * CNY) 830 : UTG = (UXG * CTX) + (UYG * CTY) 831 : RETG = CONS1.VAL(4,NLCEG) 832 : GAMG = CONS1.VAL(5,NLCEG) 833 : RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG)) 834 : RECG = RECG / ROG 835 : PG = (GAMG - 1.0D0) * (RETG - RECG) 836 : VOLG = MPOVOL.VPOCHA(NLCEG,1) 837 : C 838 : ROD = CONS1.VAL(1,NLCED) 839 : RUXD = CONS1.VAL(2,NLCED) 840 : RUYD = CONS1.VAL(3,NLCED) 841 : UXD = RUXD / ROD 842 : UYD = RUYD / ROD 843 : RETD = CONS1.VAL(4,NLCED) 844 : GAMD = CONS1.VAL(5,NLCED) 845 : RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD)) 846 : RECD = RECD / ROD 847 : PD = (GAMD - 1.0D0) * (RETD - RECD) 848 : UND = (UXD * CNX) + (UYD * CNY) 849 : UTD = (UXD * CTX) + (UYD * CTY) 850 : C 851 : IF(NLCEG .EQ. NLCED)THEN 852 : C Murs ou condition limite 853 : NLLIM=MLELIM.LECT(NGCF) 854 : IF(NLLIM .EQ.0)THEN 855 : C Mur 856 : UND = -1.0D0 * UNG 857 : UTD = UTG 858 : UXD = UND * CNX + UTD * CTX 859 : UYD = UND * CNY + UTD * CTY 860 : RUXD = ROD*UXD 861 : RUYD = ROD*UYD 862 : ELSE 863 : ROD = MPLIM.VPOCHA(NLLIM,1) 864 : UXD = MPLIM.VPOCHA(NLLIM,2) 865 : UYD = MPLIM.VPOCHA(NLLIM,3) 866 : RUXD = ROD*UXD 867 : RUYD = ROD*UYD 868 : PD = MPLIM.VPOCHA(NLLIM,4) 869 : GAMD = GAMG 870 : UND = (UXD * CNX) + (UYD * CNY) 871 : UTD = (UXD * CTX) + (UYD * CTY) 872 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD) 873 : & +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD))) 874 : ENDIF 875 : ENDIF 876 : C 877 : C**************** We check that density and pressure are positive 878 : C 879 : IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR. 880 : & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN 881 : WRITE(IOIMP,*) 'FMM' 882 : WRITE(IOIMP,*) 'ITERATION JAC',IBLJAC 883 : WRITE(IOIMP,*) 'ITERATION ICEN',ICEN 884 : WRITE(IOIMP,*) ROG,PG 885 : WRITE(IOIMP,*) ROD,PD 886 : C 887 : C 22 0 888 : C******************* Opération malvenue. Résultat douteux 889 : C 889 : CALL ERREUR(22) 890 : IPROBL = 1 891 : GOTO 9998 892 : ENDIF 893 : C 894 : C 895 : C**************** We compute the flux 896 : C FLU.VALVEC(1:4,IFAC) flux in IFAC face 897 : C 1 -> mass 898 : C 2 -> momentum (x) 899 : C 3 -> momentum (y) 900 : C 4 -> energy 901 : C 902 : RHTD = RETD+PD 903 : FR = (ROD*UND) 904 : FRUX = (ROD*UND*UXD)+(PD*CNX) 905 : FRUY = (ROD*UND*UYD)+(PD*CNY) 906 : FRET = (UND*RHTD) 907 : C 908 : COEF = SURF/(2.0D0*VOLG) 909 : C 910 : BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) 911 : & + FR * COEF 912 : BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) 913 : & + FRUX * COEF 914 : BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) 915 : & + FRUY * COEF 916 : BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) 917 : & + FRET * COEF 918 : C 919 : C**************** Viscous contribution 920 : C 921 : COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG 922 : C 923 : BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) 924 : & - (COEF * ROD) 925 : BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) 926 : & - (COEF * RUXD) 927 : BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) 928 : & - (COEF * RUYD) 929 : BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) 930 : & - (COEF * RETD) 931 : C 932 : C**************** On calcule DIS1 933 : C 934 : C DIS1 in the referred report is given by 935 : C \sum_j 936 : C \left\{ 937 : C \frac{S_j}{2 \Omega_i} ~ 938 : C \left( 939 : C \sigma_{i,j}^{0} 940 : C \mathbf{U}_{j,R}^{n+1,m,l} 941 : C \right) 942 : C \right\} 943 : C 944 : C 945 : COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG) 946 : C 947 : CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG) 948 : & + (COEF * ROD) 949 : CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG) 950 : & + (COEF * RUXD) 951 : CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG) 952 : & + (COEF * RUYD) 953 : CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG) 954 : & + (COEF * RETD) 955 : 956 : ENDDO 957 : C 958 : C************* End of the loop on the faces of the element 959 : C 960 : DO ICOMP=1,4,1 961 : D1.VALVEC(ICOMP) = 0.0D0 962 : ENDDO 963 : C 964 : COEF=0.0D0 965 : DO I1=1,4,1 966 : COEF=COEF+(Q2.VAL(I1,ICEN) 967 : & *(CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN))) 968 : ENDDO 969 : COEF=COEF*(1.0D0-PHI2V.VALVEC(ICEN)) 970 : & /PHI2V.VALVEC(ICEN) 971 : C 972 : DO I1=1,4,1 973 : D1.VALVEC(I1)=(CN1.VAL(I1,ICEN) 974 : & - CN0.VAL(I1,ICEN)) 975 : & +(COEF*Q1.VAL(I1,ICEN)) 976 : c write(6,*)D1.VALVEC(I1) 977 : ENDDO 978 : C 979 : C************* At this moment D1 = \Gamma \cdot (\Delta CN1) 980 : C 981 : DO ICOMP=1,4,1 982 : C 983 : C************* First increment 984 : C 985 : D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP) 986 : & +RES.VAL(ICOMP,ICEN) 987 : & +(BN0.VAL(ICOMP,ICEN)-BN1.VAL(ICOMP,ICEN))) 988 : & /(AN0.VALVEC(ICEN)+DN0.VALVEC(ICEN)) 989 : ENDDO 990 : C 991 : C************* At this moment, in the referred report, 992 : C D1 is (\delta \mathbf{U}')^{l+1}_i in section A3 993 : C 994 : C************* Second increment 995 : C 996 : COEF=0.0D0 997 : DO I1=1,4,1 998 : COEF=COEF+(Q2.VAL(I1,ICEN)*D1.VALVEC(I1)) 999 : ENDDO 1000 : C 1001 : C************* COEF is the scalar product between Q2 and the first 1002 : C increment. 1003 : C Let us multiply it by 1004 : C b_i^0 / (a_i^0 + b_i^0) 1005 : C 1006 : COEF=COEF*AN0.VALVEC(ICEN) 1007 : & *(PHI2V.VALVEC(ICEN)-1.0D0) 1008 : COEF=COEF/(AN0.VALVEC(ICEN)+(DN0.VALVEC(ICEN) 1009 : & *PHI2V.VALVEC(ICEN))) 1010 : C 1011 : DO I1=1,4,1 1012 : D1.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN) 1013 : ENDDO 1014 : C 1015 : C************* D1 = \delta \mathbf{U} 1016 : C 1017 : DO ICOMP=1,4,1 1018 : DU.VAL(ICOMP,ICEN) = D1.VALVEC(ICOMP) 1019 : C 1020 : CONS1.VAL(ICOMP,ICEN) = CONS0.VAL(ICOMP,ICEN) 1021 : & + D1.VALVEC(ICOMP) 1022 : C 1023 : BN1.VAL(ICOMP,ICEN) = 0.0D0 1024 : CN1.VAL(ICOMP,ICEN) = 0.0D0 1025 : ENDDO 1026 : C 1027 : C********** End of the loop on the elements 1028 : ENDDO 1029 : C 1030 : C******* End of the loop on the Domains 1031 : ENDDO 1032 : C 1033 : C**** End of the Jacobi loop 1034 : ENDDO 1035 : C 1036 : 9998 CONTINUE 1037 : DO ICEN=1,NCEN,1 1038 : DO ICOMP=1,4,1 1039 : MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN) 1040 : ENDDO 1041 : ENDDO 1042 : C 1043 : SEGSUP DU 1044 : SEGSUP RES 1045 : SEGSUP CONS0 1046 : SEGSUP CONS1 1047 : SEGSUP Q1 1048 : SEGSUP Q2 1049 : SEGSUP SIGMA0 1050 : SEGSUP AN0 1051 : SEGSUP AN0 1052 : SEGSUP BN0 1053 : SEGSUP BN1 1054 : SEGSUP CN0 1055 : SEGSUP CN1 1056 : SEGSUP DN0 1057 : SEGSUP PHI2V 1058 : SEGSUP USDTI 1059 : C 1060 : SEGDES MPOVSU 1061 : SEGDES MPOVDI 1062 : SEGDES MPOVOL 1063 : SEGDES MPNORM 1064 : SEGDES MPOCO 1065 : SEGDES MPCOEV 1066 : C 1067 : SEGSUP MPOCO1 1068 : C 1069 : SEGDES MPODU 1070 : C 1071 : SEGDES MELEFL 1072 : SEGSUP MLENT1 1073 : SEGSUP MLENT2 1074 : C 1075 : SEGSUP MLELIM 1076 : IF(MPLIM .GT. 0) SEGDES MPLIM 1077 : C 1078 : DO ISOUP=1,NSOUPO,1 1079 : IPT2=MLESUP.LECT(ISOUP) 1080 : SEGDES IPT2 1081 : ENDDO 1082 : IPT2 = MELTFA 1083 : IF(ISOUP .GT. 1) SEGDES IPT2 1084 : SEGSUP MLESUP 1085 : C 1086 : 9999 CONTINUE 1087 : RETURN 1088 : END 1089 : C 1090 : 1091 : 1092 :
© Cast3M 2003 - All rights reserved.
Disclaimer