1 : C KFM11 SOURCE BECC 06/05/15 21:15:02 5442 2 : SUBROUTINE KFM11(IRES,IRN,IGN,IRETN,IGAMN,IVCO, 3 : & ICHPSU,ICHPDI,ICHPVO,INORM, 4 : & MELEMC,MELEMF,MELEFE,DTPS,DCFL, 5 : & MELLIM,ICHLIM,NJAC,ICOEV, 6 : & IDU,IPROBL) 7 : C************************************************************************ 8 : C 9 : C PROJET : CASTEM 2000 10 : C 11 : C NOM : KFM11 12 : C 13 : C DESCRIPTION : Voir aussi KFM1 14 : C 14 : C Cas deux dimensions, gaz "calorically perfect" 15 : C Méthode de ralaxation de type Jacobi par point 16 : C 17 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 18 : C 19 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF 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 MELEFE : 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 : Avril 2002 : creation 88 : C Janvier 2003: implementation de condition aux limites 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 : CC & ,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 : & ,MELEMC,MELEMF,MELEFE,IDU,NFAC 117 : & ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED 118 : & ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP 119 : & ,I1,I2, IBLJAC, IPROBL, IVCO, ICOEV 120 : REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, REC, PG, VOLG, DIAMG 121 : & , ROD, RUXD, RUYD, UND, RETD, GAMD, PD, VOLD, DIAMD 122 : & , CNX, CNY, DCFL, DTPS 123 : & , SURF, SIGMAF 124 : & , UNSDT, UNSDTF, UXD, UYD 125 : & , FR, FRUX, FRUY, FRET, COEF 126 : & , CTX, CTY, RHTD, RHTG, UTD, UTG, UXG, UYG, QG, CG, UCOG 127 : & , ROF,PF,UNF, GAMF, UCOF, CF, PHI2, COEF1 128 : & , QN 129 : CHARACTER*8 TYPE 130 : C 131 : C**** LES INCLUDES 132 : C 133 : -INC CCOPTIO 134 : -INC SMCHPOI 135 : POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL 136 : & , MPODU.MPOVAL, MPOVOL.MPOVAL 137 : & , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL 138 : & , MPNORM.MPOVAL, MPLIM.MPOVAL, MPORES.MPOVAL, MPOCO.MPOVAL 139 : & , MPOCO1.MPOVAL, MPCOEV.MPOVAL 140 : -INC SMCOORD 141 : -INC SMELEME 142 : -INC SMLMOTS 143 : -INC SMLENTI 144 : POINTEUR MLELIM.MLENTI 145 : C 146 : C**** POINTEUR vecteurs 147 : C 148 : SEGMENT VEC 149 : REAL*8 VALVEC(I1) 150 : ENDSEGMENT 151 : POINTEUR D1.VEC,SIGMA0.VEC, AN0.VEC, DN0.VEC 152 : & ,USDTI.VEC, PHI2V.VEC, SIGMAV.VEC 153 : C & ,D2.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 233 : SEGINI, MPOCO1=MPOCO 234 : C 235 : C**** MPOCO1 is the "corrected cut-off" speed. 236 : C It is bigger than the given cut-off and the local 237 : C speed. 238 : C 239 : IF(IGEOMF .NE. MELEMF)THEN 240 : WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso' 241 : WRITE(IOIMP,*) 'Probleme de SPG' 242 : C 21 2 243 : C Données incompatibles 243 : CALL ERREUR(21) 244 : GOTO 9999 245 : ENDIF 246 : IF(IGEOMC .NE. MELEMC)THEN 247 : WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso' 248 : WRITE(IOIMP,*) 'Probleme de SPG' 249 : C 21 2 250 : C Données incompatibles 250 : CALL ERREUR(21) 251 : GOTO 9999 252 : ENDIF 253 : C 254 : C**** On initialise GAMMA=I+Q1.(Q2)^t 255 : C 256 : I1=4 257 : I2=NCEN 258 : SEGINI Q1 259 : SEGINI Q2 260 : C 261 : C**** On initialise le segment de travail D1 et D2 262 : C 263 : I1=4 264 : SEGINI D1 265 : C SEGINI D2 266 : C 267 : C**** On initialise AN0, DN0, PHI2V 268 : C 269 : I1=NCEN 270 : SEGINI AN0 271 : SEGINI DN0 272 : SEGINI PHI2V 273 : C 274 : C**** On initialise BN0 275 : C 276 : I1=4 277 : I2=NCEN 278 : SEGINI BN0 279 : C 280 : C**** On initialise BN1 281 : C 282 : I1=4 283 : I2=NCEN 284 : SEGINI BN1 285 : C 286 : C**** On initialise CN0 287 : C 288 : I1=4 289 : I2=NCEN 290 : SEGINI CN0 291 : C 292 : C**** On initialise CN1 293 : C 294 : I1=4 295 : I2=NCEN 296 : SEGINI CN1 297 : C 298 : C**** On initialise SIGMA0, SIGMAV 299 : C 300 : I1=NFAC 301 : SEGINI SIGMA0 302 : SEGINI SIGMAV 303 : C 304 : C**** RES 305 : C 306 : I1=4 307 : I2=NCEN 308 : SEGINI RES 309 : DO ICEN=1,NCEN,1 310 : DO ICOMP=1,4,1 311 : RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP) 312 : ENDDO 313 : ENDDO 314 : SEGDES MPORES 315 : C 316 : C***** DU (increment des variables conservatives) 317 : C 318 : I1=4 319 : I2=NCEN 320 : SEGINI DU 321 : C 322 : C**** Les variables conservatives 323 : C 324 : C CONS1.VAL contient les variables conservatives + gamma 325 : C i.e. rho(i) = CONS.VAL(1,i) 326 : C rhoux(i) = CONS.VAL(2,i) 327 : C rhouy(i) = CONS.VAL(3,i) 328 : C rhoet(i) = CONS.VAL(4,i) 329 : C gam(i) = CONS.VAL(5,i) 330 : C 331 : I1=5 332 : I2=NCEN 333 : SEGINI CONS1 334 : DO ICEN=1,NCEN,1 335 : CONS1.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1) 336 : CONS1.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1) 337 : CONS1.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2) 338 : CONS1.VAL(4,ICEN)=MPRETN.VPOCHA(ICEN,1) 339 : CONS1.VAL(5,ICEN)=MPGAMN.VPOCHA(ICEN,1) 340 : ENDDO 341 : SEGDES MPRN 342 : SEGDES MPRETN 343 : SEGDES MPGN 344 : SEGDES MPGAMN 345 : C 346 : SEGINI, CONS2=CONS1 347 : C 348 : C**** BOUCLE JACOBI 349 : C 350 : DO IBLJAC=1,NJAC,1 351 : C 352 : C**** BOUCLE SUR FACEL pour le calcul de AN0, BN0 353 : C 354 : DO NLCF = 1, NFAC 355 : C 356 : C******* NLCF = numero local du centre de facel 357 : C NGCF = numero global du centre de facel 358 : C NLCF1 = numero local du centre de face 359 : C NGCEG = numero global du centre ELT "gauche" 360 : C NLCEG = numero local du centre ELT "gauche" 361 : C NGCED = numero global du centre ELT "droite" 362 : C NLCED = numero local du centre ELT "droite" 363 : C 364 : NGCEG = IPT2.NUM(1,NLCF) 365 : NGCED = IPT2.NUM(3,NLCF) 366 : NGCF = IPT2.NUM(2,NLCF) 367 : NLCF1 = MLENT2.LECT(NGCF) 368 : NLCEG = MLENT1.LECT(NGCEG) 369 : NLCED = MLENT1.LECT(NGCED) 370 : C 371 : C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris. 372 : C 373 : IF(NLCF .NE. NLCF1)THEN 374 : WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso' 375 : WRITE(IOIMP,*) 'Probleme de SPG' 376 : C 21 2 377 : C Données incompatibles 377 : CALL ERREUR(21) 378 : GOTO 9999 379 : ENDIF 380 : C 381 : CNX = MPNORM.VPOCHA(NLCF,1) 382 : CNY = MPNORM.VPOCHA(NLCF,2) 383 : CTX = -1.0D0 * CNY 384 : CTY = CNX 385 : SURF= MPOVSU.VPOCHA(NLCF,1) 386 : C 387 : C******* Recuperation des Etats "gauche" et "droite" 388 : C 389 : C 390 : ROG = CONS2.VAL(1,NLCEG) 391 : RUXG = CONS2.VAL(2,NLCEG) 392 : RUYG = CONS2.VAL(3,NLCEG) 393 : UXG = RUXG / ROG 394 : UYG = RUYG / ROG 395 : UNG = (UXG * CNX) + (UYG * CNY) 396 : UTG = (UXG * CTX) + (UYG * CTY) 397 : RETG = CONS2.VAL(4,NLCEG) 398 : GAMG = CONS2.VAL(5,NLCEG) 399 : REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG)) 400 : REC = REC / ROG 401 : PG = (GAMG - 1.0D0) * (RETG - REC) 402 : VOLG = MPOVOL.VPOCHA(NLCEG,1) 403 : DIAMG = MPOVDI.VPOCHA(NLCEG,1) 404 : C 405 : ROD = CONS2.VAL(1,NLCED) 406 : RUXD = CONS2.VAL(2,NLCED) 407 : RUYD = CONS2.VAL(3,NLCED) 408 : UXD = RUXD / ROD 409 : UYD = RUYD / ROD 410 : RETD = CONS2.VAL(4,NLCED) 411 : GAMD = CONS2.VAL(5,NLCED) 412 : REC = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD)) 413 : REC = REC / ROD 414 : PD = (GAMD - 1.0D0) * (RETD - REC) 415 : VOLD = MPOVOL.VPOCHA(NLCED,1) 416 : DIAMD = MPOVDI.VPOCHA(NLCED,1) 417 : UND = (UXD * CNX) + (UYD * CNY) 418 : UTD = (UXD * CTX) + (UYD * CTY) 419 : C 420 : IF(NLCEG .EQ. NLCED)THEN 421 : C Murs ou condition limite 422 : NLLIM=MLELIM.LECT(NGCF) 423 : IF(NLLIM .EQ.0)THEN 424 : C Mur 425 : UND = -1.0D0 * UNG 426 : UTD = UTG 427 : UXD = UND * CNX + UTD * CTX 428 : UYD = UND * CNY + UTD * CTY 429 : RUXD= ROD*UXD 430 : RUYD= ROD*UYD 431 : ELSE 432 : ROD = MPLIM.VPOCHA(NLLIM,1) 433 : UXD = MPLIM.VPOCHA(NLLIM,2) 434 : UYD = MPLIM.VPOCHA(NLLIM,3) 435 : RUXD= ROD*UXD 436 : RUYD= ROD*UYD 437 : PD = MPLIM.VPOCHA(NLLIM,4) 438 : GAMD = GAMG 439 : UND = (UXD * CNX) + (UYD * CNY) 440 : UTD = (UXD * CTX) + (UYD * CTY) 441 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD)+ 442 : & (0.5D0*ROD*((UXD*UXD)+(UYD*UYD))) 443 : VOLD = VOLG 444 : ENDIF 445 : ENDIF 446 : C 447 : C********** We check that density and pressure are positive 448 : C 449 : IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR. 450 : & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN 451 : WRITE(IOIMP,*) 'FMM' 452 : C 22 0 453 : C*********** Opération malvenue. Résultat douteux 454 : C 454 : CALL ERREUR(22) 455 : IPROBL = 1 456 : GOTO 9998 457 : ENDIF 458 : C 459 : RHTG=RETG+PG 460 : RHTD=RETD+PD 461 : C 462 : C********** FIRST JACOBI ITERATION 463 : C 464 : IF(IBLJAC .EQ. 1)THEN 465 : C 466 : C************* We compute the interfacial state 467 : C 468 : ROF=0.5D0*(ROG+ROD) 469 : PF=0.5D0*(PG+PD) 470 : UNF=0.5D0*(UNG+UND) 471 : GAMF=0.5D0*(GAMG+GAMD) 472 : C 473 : C************* We compute SIGMA and the corrected cut-off speed 474 : C Here SIGMA is the spectral of the preconditioned Roe 475 : C matrix (|\Gamma^{-1} A| in the referred report) 476 : C 477 : UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1)) 478 : UCOF=MAX(UCOF,(((UNG*UNG)+(UTG*UTG))**0.5D0)) 479 : UCOF=MAX(UCOF,(((UND*UND)+(UTD*UTD))**0.5D0)) 480 : IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1)) MPOCO1.VPOCHA(NLCEG 481 : $ ,1)= UCOF 482 : IF(UCOF .GT. MPOCO1.VPOCHA(NLCED,1)) MPOCO1.VPOCHA(NLCED 483 : $ ,1)= UCOF 484 : C 485 : QN=ABS(UNF) 486 : CF=(GAMF*PF/ROF)**0.5D0 487 : IF(UCOF .GE. CF)THEN 488 : PHI2=1.0D0 489 : ELSEIF(QN .GT. CF)THEN 490 : PHI2=1.0D0 491 : ELSEIF(QN .LT. UCOF)THEN 492 : PHI2=UCOF / CF 493 : ELSE 494 : PHI2 = QN / CF 495 : ENDIF 496 : PHI2=PHI2*PHI2 497 : C 498 : SIGMAF=(1.0D0-PHI2)*QN 499 : SIGMAF=SIGMAF*SIGMAF 500 : SIGMAF=SIGMAF+(4.0D0*PHI2*CF*CF) 501 : SIGMAF=SIGMAF**0.5D0 502 : SIGMAF=SIGMAF+((1.0D0+PHI2)*QN) 503 : SIGMAF=0.5D0*SIGMAF 504 : c 505 : c write(*,*) qn,sigmaf,phi2 506 : c 507 : c write(*,*) 'UCOF',UCOF 508 : c write(*,*) 'PHI2',PHI2 509 : c write(*,*) 'JACOBI,SIGMAF',IBLJAC,SIGMAF 510 : c 511 : C 512 : C************* We store SIGMA into SIGMA0 513 : C 514 : SIGMA0.VALVEC(NLCF)=SIGMAF 515 : C 516 : C************* We compute the dual time step 517 : C 518 : UNSDT=USDTI.VALVEC(NLCEG) 519 : UNSDTF = SIGMAF / DIAMG 520 : IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG)=UNSDTF 521 : UNSDT=USDTI.VALVEC(NLCED) 522 : UNSDTF = SIGMAF / DIAMD 523 : IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCED)=UNSDTF 524 : C 525 : C************* We compute AN0 526 : C ANO in the referred report is equal to 527 : C \left( 528 : C \frac{1}{\Delta \tau_i ^0} 529 : C + \frac{1}{2{\Omega_i}} 530 : C \sum_j \sigma_{i,j}^{0} S_j 531 : C \right) 532 : C 533 : AN0.VALVEC(NLCEG)=AN0.VALVEC(NLCEG) + 534 : & (0.5D0*SIGMAF*SURF/VOLG) 535 : IF(NLCEG .NE. NLCED) 536 : & AN0.VALVEC(NLCED)=AN0.VALVEC(NLCED) + 537 : & (0.5D0*SIGMAF*SURF/VOLD) 538 : C 539 : C************* We compute SIGMAV 540 : C 541 : C COEF=FG 542 : COEF=(MCOORD.XCOOR(((NGCEG-1)*3)+1)- 543 : & MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX 544 : COEF=COEF+ 545 : & ((MCOORD.XCOOR(((NGCEG-1)*3)+2)- 546 : & MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY) 547 : C COEF1=FD 548 : COEF1=(MCOORD.XCOOR(((NGCED-1)*3)+1)- 549 : & MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX 550 : COEF1=COEF1+ 551 : & ((MCOORD.XCOOR(((NGCED-1)*3)+2)- 552 : & MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY) 553 : C COEF=|FG|+|FD| 554 : COEF=ABS(COEF)+ABS(COEF1) 555 : SIGMAF=0.5D0*(MPCOEV.VPOCHA(NLCEG,1)+ 556 : & MPCOEV.VPOCHA(NLCED,1))/COEF 557 : SIGMAV.VALVEC(NLCF)=SIGMAF 558 : C 559 : C write(*,*) 'sigmaf =',sigmaf 560 : C 561 : C************* We compute the contribution of SIGMAV to 562 : C the diagonal part of the matrix to invert 563 : C 564 : C In the given report, DNO is a_i^0, i.e. 565 : C DN0 = AN0 (previously defined) + 566 : C + (\frac{1}{{\Omega_i}} 567 : C \sum_j \sigma_{V,i,j}^{0} S_j) 568 : C + 1 / \Delta t^n 569 : C 570 : COEF=SURF/VOLG 571 : DN0.VALVEC(NLCEG)=DN0.VALVEC(NLCEG)+COEF*SIGMAF 572 : IF(NLCEG .NE. NLCED)THEN 573 : COEF=SURF/VOLD 574 : DN0.VALVEC(NLCED)=DN0.VALVEC(NLCED)+COEF*SIGMAF 575 : ENDIF 576 : ENDIF 577 : C 578 : C********** ALL JACOBI ITERATIONS 579 : C 580 : C 581 : C********** On calcule BN1 582 : C 583 : C BN1 in the referred report is given by 584 : C \sum_j 585 : C \left\{ 586 : C \frac{S_j}{2 \Omega_i} ~ 587 : C \left( 588 : C \mathcal{F}(\mathbf{U}_{j,R}^{\ell}) 589 : C - 590 : C 2 \sigma_{V,i,j}^{0} 591 : C \mathbf{U}_{j,R}^{\ell} 592 : C _right) 593 : C \right\} 594 : C 595 : C Central contribution 596 : C 597 : FR=(ROD*UND) 598 : FRUX=(ROD*UND*UXD)+(PD*CNX) 599 : FRUY=(ROD*UND*UYD)+(PD*CNY) 600 : FRET=(UND*RHTD) 601 : C 602 : COEF=SURF/(2.0D0*VOLG) 603 : BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) + 604 : & (COEF * FR) 605 : BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) + 606 : & (COEF * FRUX) 607 : BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) + 608 : & (COEF * FRUY) 609 : BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) + 610 : & (COEF * FRET) 611 : C 612 : C Viscous contribution 613 : C 614 : COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG 615 : BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) - 616 : & (COEF * ROD) 617 : BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) - 618 : & (COEF * RUXD) 619 : BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) - 620 : & (COEF * RUYD) 621 : BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) - 622 : & (COEF * RETD) 623 : C 624 : C********** On calcule CN1 625 : C 626 : C CN1 in the referred report is given by 627 : C \sum_j 628 : C \left\{ 629 : C \frac{S_j}{2 \Omega_i} ~ 630 : C \left( 631 : C \sigma_{i,j}^{0} 632 : C \mathbf{U}_{j,R}^{\ell} 633 : C \right) 634 : C \right\} 635 : C 636 : C 637 : COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG) 638 : CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG) + 639 : & (COEF * ROD) 640 : CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG) + 641 : & (COEF * RUXD) 642 : CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG) + 643 : & (COEF * RUYD) 644 : CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG) + 645 : & (COEF * RETD) 646 : C 647 : IF(NLCEG .NE. NLCED)THEN 648 : C 649 : C********** Domain interieur 650 : C 651 : C Central contribution 652 : C 653 : FR=(ROG*UNG) 654 : FRUX=(ROG*UNG*UXG)+(PG*CNX) 655 : FRUY=(ROG*UNG*UYG)+(PG*CNY) 656 : FRET=(UNG*RHTG) 657 : COEF=SURF/(2.0D0*VOLD) 658 : C 659 : BN1.VAL(1,NLCED) = BN1.VAL(1,NLCED) - 660 : & (COEF * FR) 661 : BN1.VAL(2,NLCED) = BN1.VAL(2,NLCED) - 662 : & (COEF * FRUX) 663 : BN1.VAL(3,NLCED) = BN1.VAL(3,NLCED) - 664 : & (COEF * FRUY) 665 : BN1.VAL(4,NLCED) = BN1.VAL(4,NLCED) - 666 : & (COEF * FRET) 667 : C 668 : C Viscous contribution 669 : C 670 : COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLD 671 : BN1.VAL(1,NLCED) = BN1.VAL(1,NLCED) - 672 : & (COEF * ROG) 673 : BN1.VAL(2,NLCED) = BN1.VAL(2,NLCED) - 674 : & (COEF * RUXG) 675 : BN1.VAL(3,NLCED) = BN1.VAL(3,NLCED) - 676 : & (COEF * RUYG) 677 : BN1.VAL(4,NLCED) = BN1.VAL(4,NLCED) - 678 : & (COEF * RETG) 679 : C 680 : C Numerical diffusity contribution 681 : C 682 : COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLD) 683 : CN1.VAL(1,NLCED) = CN1.VAL(1,NLCED) + 684 : & (COEF * ROG) 685 : CN1.VAL(2,NLCED) = CN1.VAL(2,NLCED) + 686 : & (COEF * RUXG) 687 : CN1.VAL(3,NLCED) = CN1.VAL(3,NLCED) + 688 : & (COEF * RUYG) 689 : CN1.VAL(4,NLCED) = CN1.VAL(4,NLCED) + 690 : & (COEF * RETG) 691 : 692 : ENDIF 693 : C 694 : C********** Fin boucle sur les faces 695 : C 696 : ENDDO 697 : 698 : C 699 : C********** Loop on the centre 700 : C Computation of GAMMA 701 : C DMAT 702 : C 703 : DO ICEN=1,NCEN,1 704 : C 705 : IF(IBLJAC .EQ. 1)THEN 706 : C 707 : ROG = CONS2.VAL(1,ICEN) 708 : RUXG = CONS2.VAL(2,ICEN) 709 : RUYG = CONS2.VAL(3,ICEN) 710 : UXG=RUXG/ROG 711 : UYG=RUYG/ROG 712 : RETG = CONS2.VAL(4,ICEN) 713 : GAMG = CONS2.VAL(5,ICEN) 714 : REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG)) 715 : REC = REC / ROG 716 : PG = (GAMG - 1.0D0) * (RETG - REC) 717 : C 718 : C************* Note that the Cut-off is the corrected one 719 : C 720 : CG=(GAMG*PG/ROG)**0.5D0 721 : UCOG=MPOCO1.VPOCHA(ICEN,1) 722 : QG=2.0D0*REC/ROG 723 : QG=QG**0.5D0 724 : C 725 : IF(UCOG .GE. CG)THEN 726 : PHI2=1.0D0 727 : ELSEIF(QG .GT. CG)THEN 728 : PHI2=1.0D0 729 : ELSEIF(QG .LT. UCOG)THEN 730 : PHI2=UCOG / CG 731 : ELSE 732 : PHI2 = QG / CG 733 : ENDIF 734 : PHI2=PHI2*PHI2 735 : PHI2V.VALVEC(ICEN)=PHI2 736 : C 737 : C************* We compute GAMMA 738 : C GAMMA=I+((1-PHI2)/PHI2*Q1.(Q2)^t) 739 : C 740 : C In the referred report, Q2 is the line vector 741 : C 742 : C ( \frac{q^2}{2} , -u , -v , 1) 743 : C 744 : C and Q1 is the column vector 745 : C 746 : C \frac{\gamma - 1}{c^2} 747 : C \left( 748 : C \begin{array}{c} 749 : C 1 \\ 750 : C u \\ 751 : C v \\ 752 : C H \\ 753 : C \end{array} 754 : C 755 : C 756 : Q2.VAL(1,ICEN)=(0.5D0*QG*QG) 757 : Q2.VAL(2,ICEN)=-1.0D0*RUXG/ROG 758 : Q2.VAL(3,ICEN)=-1.0D0*RUYG/ROG 759 : Q2.VAL(4,ICEN)=1.0D0 760 : C 761 : COEF1=(GAMG-1.0D0)/(CG*CG) 762 : COEF=(1.0D0/COEF1)+(0.5D0*QG*QG) 763 : C COEF=HT 764 : Q1.VAL(1,ICEN)=COEF1 765 : Q1.VAL(2,ICEN)=COEF1*RUXG/ROG 766 : Q1.VAL(3,ICEN)=COEF1*RUYG/ROG 767 : Q1.VAL(4,ICEN)=COEF1*COEF 768 : C 769 : AN0.VALVEC(ICEN)=USDTI.VALVEC(ICEN)*2.0D0/DCFL 770 : & +AN0.VALVEC(ICEN) 771 : DN0.VALVEC(ICEN)=DN0.VALVEC(ICEN)+(1.0D0/DTPS) 772 : C 773 : C 774 : C************* We save BN1, CN1 into BN0 and CN0 775 : C 776 : DO ICOMP=1,4,1 777 : BN0.VAL(ICOMP,ICEN) = BN1.VAL(ICOMP,ICEN) 778 : ENDDO 779 : C 780 : DO ICOMP=1,4,1 781 : CN0.VAL(ICOMP,ICEN) = CN1.VAL(ICOMP,ICEN) 782 : ENDDO 783 : C 784 : ENDIF 785 : C 786 : C********** IBLJAC => 1 787 : C 788 : DO ICOMP=1,4,1 789 : D1.VALVEC(ICOMP) = 0.0D0 790 : C D2.VALVEC(ICOMP) = 0.0D0 791 : ENDDO 792 : C 793 : COEF=0.0D0 794 : DO I1=1,4,1 795 : COEF=COEF+(Q2.VAL(I1,ICEN)*(CN1.VAL(I1,ICEN) - CN0.VAL(I1 796 : $ ,ICEN))) 797 : ENDDO 798 : COEF=COEF*(1.0D0-PHI2V.VALVEC(ICEN))/PHI2V.VALVEC(ICEN) 799 : C 800 : DO I1=1,4,1 801 : D1.VALVEC(i1)=(CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN))+ 802 : & (COEF*Q1.VAL(I1,ICEN)) 803 : ENDDO 804 : C 805 : C********** At this moment D1 = \Gamma \cdot (\Delta CN1) 806 : C 807 : DO ICOMP=1,4,1 808 : C 809 : C************* First increment 810 : C 811 : D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP)+ 812 : & RES.VAL(ICOMP,ICEN)+(BN0.VAL(ICOMP,ICEN) - 813 : & BN1.VAL(ICOMP,ICEN)))/ 814 : & (AN0.VALVEC(ICEN)+DN0.VALVEC(ICEN)) 815 : ENDDO 816 : C 817 : C********** At this moment, in the referred report, 818 : C D1 is (\delta \mathbf{U}')^{l+1}_i in section A3 819 : C 820 : C********** Second increment 821 : C 822 : COEF=0.0D0 823 : DO I1=1,4,1 824 : COEF=COEF+(Q2.VAL(I1,ICEN)*D1.VALVEC(I1)) 825 : ENDDO 826 : C 827 : C********** COEF is the scalar product between Q2 and the first 828 : C increment. 829 : C Let us multiply it by 830 : C b_i^0 / (a_i^0 + b_i^0) 831 : C 832 : COEF=COEF*AN0.VALVEC(ICEN)*(PHI2V.VALVEC(ICEN)-1.0D0) 833 : COEF=COEF/(AN0.VALVEC(ICEN)+(DN0.VALVEC(ICEN) 834 : & *PHI2V.VALVEC(ICEN))) 835 : C 836 : DO I1=1,4,1 837 : C D2.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN) 838 : D1.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN) 839 : ENDDO 840 : C 841 : C********** D2 = \delta \mathbf{U} 842 : C 843 : DO ICOMP=1,4,1 844 : C DU.VAL(ICOMP,ICEN)=D2.VALVEC(ICOMP) 845 : C CONS2.VAL(ICOMP,ICEN) = CONS1.VAL(ICOMP,ICEN) + 846 : C & D2.VALVEC(ICOMP) 847 : DU.VAL(ICOMP,ICEN)=D1.VALVEC(ICOMP) 848 : CONS2.VAL(ICOMP,ICEN) = CONS1.VAL(ICOMP,ICEN) + 849 : & D1.VALVEC(ICOMP) 850 : BN1.VAL(ICOMP,ICEN)=0.0D0 851 : CN1.VAL(ICOMP,ICEN)=0.0D0 852 : ENDDO 853 : ENDDO 854 : C 855 : C**** Fin boucle JACOBI 856 : C 857 : ENDDO 858 : C 859 : 9998 CONTINUE 860 : c IF(IPROBL .EQ. 0)THEN 861 : DO ICEN=1,NCEN,1 862 : DO ICOMP=1,4,1 863 : MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN) 864 : ENDDO 865 : ENDDO 866 : c ELSE 867 : c DO ICEN=1,NCEN,1 868 : c DO ICOMP=1,4,1 869 : c MPODU.VPOCHA(ICEN,ICOMP) = 0.0D0 870 : c ENDDO 871 : c ENDDO 872 : c ENDIF 873 : C 874 : SEGSUP DU 875 : SEGSUP RES 876 : SEGSUP CONS1 877 : SEGSUP CONS2 878 : SEGSUP Q1 879 : SEGSUP Q2 880 : SEGSUP D1 881 : C SEGSUP D2 882 : SEGSUP SIGMA0 883 : SEGSUP SIGMAV 884 : SEGSUP AN0 885 : SEGSUP BN0 886 : SEGSUP BN1 887 : SEGSUP CN0 888 : SEGSUP CN1 889 : SEGSUP DN0 890 : SEGSUP PHI2V 891 : SEGSUP USDTI 892 : C 893 : SEGDES MPOVSU 894 : SEGDES MPOVDI 895 : SEGDES MPOVOL 896 : SEGDES MPNORM 897 : SEGDES MPOCO 898 : SEGDES MPCOEV 899 : C 900 : SEGSUP MPOCO1 901 : C 902 : SEGDES MPODU 903 : C 904 : SEGDES IPT2 905 : SEGSUP MLENT1 906 : SEGSUP MLENT2 907 : C 908 : SEGSUP MLELIM 909 : IF(MPLIM .GT. 0) SEGDES MPLIM 910 : C 911 : 9999 CONTINUE 912 : RETURN 913 : END 914 : C 915 : 916 : 917 :
© Cast3M 2003 - All rights reserved.
Disclaimer