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