1 : C KFM14 SOURCE BECC 06/05/15 21:15:28 5442 2 : C KFM14 SOURCE CHAT 05/01/13 00:55:31 5004 3 : SUBROUTINE KFM14(IRES,IRN,IGN,IRETN,IGAMN,IVCO, 4 : & ICHPSU,ICHPDI,ICHPVO,INORM, 5 : & MELEMC,MELEMF,MELEFL,MELTFA,DTPS,DCFL, 6 : & MELLIM,ICHLIM,NJAC,ICOEF,ICOEV, 7 : & IDU,IPROBL) 8 : C************************************************************************ 9 : C 10 : C PROJET : CASTEM 2000 11 : C 12 : C NOM : KFM14 13 : C 14 : C DESCRIPTION : Voir aussi KFM1 et KFM12 15 : C 15 : C Cas 3D, gaz "calorically perfect" 16 : C relaxation de type SGS 17 : C 18 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 19 : C 20 : C AUTEUR : T. KLOCZKO, DEN/DM2S/SFME/LTMF 21 : C 22 : C************************************************************************ 23 : C ENTREES 24 : C 25 : C 26 : C 1) Pointeurs des CHPOINTs 27 : C 28 : C IRES : CHPOINT 'CENTRE' contenant le residu 29 : C 30 : C IRN : CHPOINT 'CENTRE' contenant la masse volumique 31 : C 32 : C IGN : CHPOINT 'CENTRE' contenant la q.d.m. 33 : C 34 : C IRETN : CHPOINT 'CENTRE' contenant l'energie totale 35 : C 36 : C IGAMN : CHPOINT 'CENTRE' contenant le gamma 37 : C 38 : C IVCO : CHPOINT 'CENTRE' contenant le cut-off 39 : C 40 : C ICHLIM : CHPOINT conditions aux bords 41 : C 42 : C ICOEV : CHPOINT 'CENTRE' contenant le (gamma - 1) K / (R \rho) 43 : C coeff pour calculer le ray. spect. de la matrice 44 : C visc. 45 : C 46 : C 3) Pointeurs de CHPOINTs de la table DOMAINE 47 : C 48 : C ICHPSU : CHPOINT "FACE" contenant la surface des faces 49 : C 50 : C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum 51 : C de chaque element 52 : C 53 : C ICHPVO : CHPOINT "CENTRE" contenant le volume 54 : C de chaque element 55 : C 56 : C INORM : CHPOINT "CENTRE" contenant le normales aux faces 57 : C 58 : C 59 : C 4) Pointeurs de MELEME de la table DOMAINE 60 : C 61 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES 62 : C 63 : C MELEMF : MELEME 'FACE' du SPG des FACES 64 : C 65 : C MELEFL : MELEME 'FACEL' du connectivité FACES -> ELEM 66 : C 67 : C MELLIM : MELEME SPG conditions aux bords 68 : C 69 : C 5) 70 : C 71 : C DCFL = le double de la CFL 72 : C 73 : C DTPS = pas de temps physique 74 : C 75 : C NJAC = nombre d'iterations dans le jacobien 76 : C 77 : C 78 : C SORTIES 79 : C 80 : C IDU : resultat 81 : C 82 : C IPROBL : si > 0, il y a un probleme 83 : C 84 : C************************************************************************ 85 : C 86 : C HISTORIQUE (Anomalies et modifications éventuelles) 87 : C 88 : C HISTORIQUE : 2004 création 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, RUZG, RETG, GAMG, RHTG, RECG 124 : & , ROD, RUXD, RUYD, RUZD, RETD, GAMD, RHTD, RECD 125 : & , UXG, UYG, UZG, UNG, UTG, UVG, PG, VOLG, DIAMG, SURF 126 : & , UXD, UYD, UZD, UND, UTD, UVD, PD 127 : & , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ 128 : & , DCFL, DTPS, UNSDT, UNSDTF 129 : & , FR, FRUX, FRUY, FRUZ, FRET 130 : & , QG, CG, UCOG, ROF, PF, UNF, GAMF, UCOF, CF, SIGMAF 131 : & , QN, PHI2, COEF, COEF1 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 dans kfm14.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 dans kfm14.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=5 262 : I2=NCEN 263 : SEGINI Q1 264 : SEGINI Q2 265 : C 266 : C**** On initialise le segment de travail D1 267 : C 268 : I1=5 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=5 281 : I2=NCEN 282 : SEGINI BN0 283 : C 284 : C**** On initialise BN1 285 : C 286 : I1=5 287 : I2=NCEN 288 : SEGINI BN1 289 : C 290 : C**** On initialise CN0 291 : C 292 : I1=5 293 : I2=NCEN 294 : SEGINI CN0 295 : C 296 : C**** On initialise CN1 297 : C 298 : I1=5 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=5 311 : I2=NCEN 312 : SEGINI RES 313 : DO ICEN=1,NCEN,1 314 : DO ICOMP=1,5,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=5 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 rhouz(i) = CONS.VAL(4,i) 333 : C rhoet(i) = CONS.VAL(5,i) 334 : C gam(i) = CONS.VAL(6,i) 335 : C 336 : I1=6 337 : I2=NCEN 338 : SEGINI CONS0 339 : DO ICEN=1,NCEN,1 340 : CONS0.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1) 341 : CONS0.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1) 342 : CONS0.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2) 343 : CONS0.VAL(4,ICEN)=MPGN.VPOCHA(ICEN,3) 344 : CONS0.VAL(5,ICEN)=MPRETN.VPOCHA(ICEN,1) 345 : CONS0.VAL(6,ICEN)=MPGAMN.VPOCHA(ICEN,1) 346 : ENDDO 347 : SEGDES MPRN 348 : SEGDES MPRETN 349 : SEGDES MPGN 350 : SEGDES MPGAMN 351 : C 352 : SEGINI, CONS1=CONS0 353 : C 354 : C**** Storage of l-1 values of F at the l Jacobi iteration 355 : C 356 : IPT1=MELTFA 357 : SEGACT IPT1 358 : NSOUPO= IPT1.LISOUS(/1) 359 : NSOUPO=MAX(NSOUPO,1) 360 : JG=NSOUPO 361 : SEGINI MLESUP 362 : IFACT=0 363 : IF(NSOUPO .EQ. 1)THEN 364 : MLESUP.LECT(1)=IPT1 365 : NFACE=IPT1.NUM(/1) 366 : NELT=IPT1.NUM(/2) 367 : IFACT=IFACT+(NFACE*NELT) 368 : ELSE 369 : DO ISOUP = 1, NSOUPO, 1 370 : IPT2=IPT1.LISOUS(ISOUP) 371 : MLESUP.LECT(ISOUP)=IPT2 372 : SEGACT IPT2 373 : NFACE=IPT2.NUM(/1) 374 : NELT=IPT2.NUM(/2) 375 : IFACT=IFACT+(NFACE*NELT) 376 : ENDDO 377 : ENDIF 378 : C 379 : C*********************************************************************** 380 : C**** FIRST JACOBI ITERATION ******************************************* 381 : C*********************************************************************** 382 : C 383 : C 384 : C**** LOOP on ELTFA to compute: 385 : C 386 : C AN0_i = 1/(2 \Omega_i) \sum_j \sigma_i,j S_j 387 : C SIGMA0_j 388 : C The cut-off in each cell 389 : C The reciprocal of the dual time step USDTI 390 : C 391 : ICEN=0 392 : IFACG=0 393 : DO ISOUP=1,NSOUPO,1 394 : IPT2=MLESUP.LECT(ISOUP) 395 : NFACE=IPT2.NUM(/1) 396 : NELT=IPT2.NUM(/2) 397 : DO IELT=1,NELT,1 398 : ICEN=ICEN+1 399 : DO IFAC=1,NFACE,1 400 : C IFACG increase with IFAC,IELT,ISOUP 401 : IFACG = IFACG+1 402 : NGCF = IPT2.NUM(IFAC,IELT) 403 : NLCF = MLENT2.LECT(NGCF) 404 : C 405 : C************* NLCF = numero local du centre de facel 406 : C NGCF = numero global du centre de facel 407 : C NGCEG = numero global du centre ELT "gauche" 408 : C NLCEG = numero local du centre ELT "gauche" 409 : C NGCED = numero global du centre ELT "droite" 410 : C NLCED = numero local du centre ELT "droite" 411 : C 412 : NGCEG = MELEFL.NUM(1,NLCF) 413 : NGCED = MELEFL.NUM(3,NLCF) 414 : NLCEG = MLENT1.LECT(NGCEG) 415 : NLCED = MLENT1.LECT(NGCED) 416 : C 417 : CNX = MPNORM.VPOCHA(NLCF,7) 418 : CNY = MPNORM.VPOCHA(NLCF,8) 419 : CNZ = MPNORM.VPOCHA(NLCF,9) 420 : CTX = MPNORM.VPOCHA(NLCF,1) 421 : CTY = MPNORM.VPOCHA(NLCF,2) 422 : CTZ = MPNORM.VPOCHA(NLCF,3) 423 : CVX = MPNORM.VPOCHA(NLCF,4) 424 : CVY = MPNORM.VPOCHA(NLCF,5) 425 : CVZ = MPNORM.VPOCHA(NLCF,6) 426 : C 427 : C************* Left or right? 428 : C 429 : IF(NLCEG .EQ. ICEN)THEN 430 : LOGINV=.FALSE. 431 : ELSEIF(NLCED .EQ. ICEN)THEN 432 : LOGINV=.TRUE. 433 : ELSE 434 : WRITE(IOIMP,*) 'Anomalie dans kfm13.eso' 435 : WRITE(IOIMP,*) 'Probleme de SPG' 436 : C 21 2 437 : C Données incompatibles 437 : CALL ERREUR(21) 438 : GOTO 9999 439 : ENDIF 440 : C************* If right, inversion 441 : IF(LOGINV) THEN 442 : NLCED=NLCEG 443 : NLCEG=ICEN 444 : CNX = -1.0D0*CNX 445 : CNY = -1.0D0*CNY 446 : CNZ = -1.0D0*CNZ 447 : CTX = -1.0D0*CTX 448 : CTY = -1.0D0*CTY 449 : CTZ = -1.0D0*CTZ 450 : CVX = -1.0D0*CVX 451 : CVY = -1.0D0*CVY 452 : CVZ = -1.0D0*CVZ 453 : ENDIF 454 : C 455 : SURF = MPOVSU.VPOCHA(NLCF,1) 456 : C 457 : C************* Recuperation des Etats "gauche" et "droite" 458 : C 459 : C 460 : ROG = CONS0.VAL(1,NLCEG) 461 : RUXG = CONS0.VAL(2,NLCEG) 462 : RUYG = CONS0.VAL(3,NLCEG) 463 : RUZG = CONS0.VAL(4,NLCEG) 464 : RETG = CONS0.VAL(5,NLCEG) 465 : GAMG = CONS0.VAL(6,NLCEG) 466 : C 467 : UXG = RUXG / ROG 468 : UYG = RUYG / ROG 469 : UZG = RUZG / ROG 470 : C 471 : UNG = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ) 472 : UTG = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ) 473 : UVG = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ) 474 : C 475 : RECG = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2) 476 : RECG = RECG / ROG 477 : PG = (GAMG - 1.0D0) * (RETG - RECG) 478 : C 479 : VOLG = MPOVOL.VPOCHA(NLCEG,1) 480 : DIAMG = MPOVDI.VPOCHA(NLCEG,1) 481 : C 482 : ROD = CONS0.VAL(1,NLCED) 483 : RUXD = CONS0.VAL(2,NLCED) 484 : RUYD = CONS0.VAL(3,NLCED) 485 : RUZD = CONS0.VAL(4,NLCED) 486 : RETD = CONS0.VAL(5,NLCED) 487 : GAMD = CONS0.VAL(6,NLCED) 488 : C 489 : UXD = RUXD / ROD 490 : UYD = RUYD / ROD 491 : UZD = RUZD / ROD 492 : C 493 : UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ) 494 : UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ) 495 : UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ) 496 : C 497 : RECD = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2) 498 : RECD = RECD / ROD 499 : PD = (GAMD - 1.0D0) * (RETD - RECD) 500 : C 501 : IF(NLCEG .EQ. NLCED)THEN 502 : C Murs ou condition limite 503 : NLLIM=MLELIM.LECT(NGCF) 504 : IF(NLLIM .EQ.0)THEN 505 : C Mur 506 : UND = -1.0D0 * UNG 507 : UTD = UTG 508 : UVD = UVG 509 : UXD = UND * CNX + UTD * CTX + UVD * CVX 510 : UYD = UND * CNY + UTD * CTY + UVD * CVY 511 : UZD = UND * CNZ + UTD * CTZ + UVD * CVZ 512 : C 513 : RUXD = ROD * UXD 514 : RUYD = ROD * UYD 515 : RUZD = ROD * UZD 516 : C 517 : ELSE 518 : ROD = MPLIM.VPOCHA(NLLIM,1) 519 : UXD = MPLIM.VPOCHA(NLLIM,2) 520 : UYD = MPLIM.VPOCHA(NLLIM,3) 521 : UZD = MPLIM.VPOCHA(NLLIM,4) 522 : PD = MPLIM.VPOCHA(NLLIM,5) 523 : C 524 : RUXD = ROD*UXD 525 : RUYD = ROD*UYD 526 : RUZD = ROD * UZD 527 : C 528 : GAMD = GAMG 529 : C 530 : UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ) 531 : UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ) 532 : UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ) 533 : C 534 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD) 535 : & +(0.5D0*ROD*(UXD**2+UYD**2+UZD**2)) 536 : ENDIF 537 : ENDIF 538 : C 539 : C************* We check that density and pressure are positive 540 : C 541 : IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR. 542 : & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN 543 : WRITE(IOIMP,*) 'FMM' 544 : WRITE(IOIMP,*) 'TEST' 545 : C 22 0 546 : C**************** Opération malvenue. Résultat douteux 547 : C 547 : CALL ERREUR(22) 548 : IPROBL = 1 549 : GOTO 9998 550 : ENDIF 551 : C 552 : C************* We compute the interfacial state 553 : C 554 : ROF = 0.5D0*(ROG+ROD) 555 : PF = 0.5D0*(PG+PD) 556 : UNF = 0.5D0*(UNG+UND) 557 : GAMF = 0.5D0*(GAMG+GAMD) 558 : C 559 : C************* We compute SIGMA and the corrected cut-off speed 560 : C 561 : UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1)) 562 : UCOF=MAX(UCOF,((UNG**2+UTG**2+UVG**2)**0.5D0)) 563 : UCOF=MAX(UCOF,((UND**2+UTD**2+UVD**2)**0.5D0)) 564 : IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1)) 565 : & MPOCO1.VPOCHA(NLCEG,1)= UCOF 566 : C 567 : QN = ABS(UNF) 568 : CF = (GAMF * PF / ROF)**0.5D0 569 : IF(UCOF .GE. CF)THEN 570 : PHI2 = 1.0D0 571 : ELSEIF(QN .GT. CF)THEN 572 : PHI2 = 1.0D0 573 : ELSEIF(QN .LT. UCOF)THEN 574 : PHI2 = UCOF / CF 575 : ELSE 576 : PHI2 = QN / CF 577 : ENDIF 578 : PHI2 = PHI2 * PHI2 579 : C 580 : SIGMAF = (1.0D0 - PHI2) * QN 581 : SIGMAF = SIGMAF * SIGMAF 582 : SIGMAF = SIGMAF + (4.0D0 * PHI2 * CF * CF) 583 : SIGMAF = SIGMAF**0.5D0 584 : SIGMAF = SIGMAF + ((1.0D0 + PHI2) * QN) 585 : SIGMAF = 0.5D0 * SIGMAF 586 : C 587 : SIGMA0.VALVEC(NLCF) = SIGMAF 588 : C 589 : C************* We compute AN0 590 : C 591 : AN0.VALVEC(NLCEG) = AN0.VALVEC(NLCEG) 592 : & + (0.5D0 * SIGMAF * SURF / VOLG) 593 : C 594 : C************* We compute the dual time step 595 : C 596 : UNSDT = USDTI.VALVEC(NLCEG) 597 : UNSDTF = SIGMAF / DIAMG 598 : IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG) = UNSDTF 599 : C 600 : C************* We compute SIGMAV 601 : C 602 : C COEF=FG 603 : COEF = (MCOORD.XCOOR(((NGCEG-1)*4)+1) 604 : & -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX 605 : COEF = COEF 606 : & + (MCOORD.XCOOR(((NGCEG-1)*4)+2) 607 : & -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY 608 : COEF = COEF 609 : & + (MCOORD.XCOOR(((NGCEG-1)*4)+3) 610 : & -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ 611 : C 612 : C COEF1=FD 613 : COEF1 = (MCOORD.XCOOR(((NGCED-1)*4)+1) 614 : & -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX 615 : COEF1 = COEF1 616 : & + (MCOORD.XCOOR(((NGCED-1)*4)+2) 617 : & -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY 618 : COEF1 = COEF1 619 : & + (MCOORD.XCOOR(((NGCED-1)*4)+3) 620 : & -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ 621 : C 622 : C COEF=|FG|+|FD| 623 : COEF = ABS(COEF)+ABS(COEF1) 624 : SIGMAF = 0.5D0 * (MPCOEV.VPOCHA(NLCEG,1) 625 : & +MPCOEV.VPOCHA(NLCED,1)) / COEF 626 : C 627 : C 628 : SIGMAV.VALVEC(NLCF) = SIGMAF 629 : C 630 : C************* We compute the contribution of SIGMAV to 631 : C the diagonal part of the matrix to invert 632 : C 633 : C In the given report, DNO is a_i^0, i.e. 634 : C DN0 = AN0 (previously defined) + 635 : C + (\frac{1}{{\Omega_i}} 636 : C \sum_j \sigma_{V,i,j}^{0} S_j) 637 : C + 1 / \Delta t^n 638 : C 639 : COEF = SURF / VOLG 640 : DN0.VALVEC(NLCEG) = DN0.VALVEC(NLCEG) + COEF * SIGMAF 641 : C 642 : C 643 : C************* We compute the flux 644 : C FLU.VALVEC(1:4,IFAC) flux in IFAC face 645 : C 1 -> mass 646 : C 2 -> momentum (x) 647 : C 3 -> momentum (y) 648 : C 4 -> energy 649 : C 650 : RHTD = RETD + PD 651 : FR = (ROD * UND) 652 : FRUX = (ROD * UND * UXD) + (PD * CNX) 653 : FRUY = (ROD * UND * UYD) + (PD * CNY) 654 : FRUZ = (ROD * UND * UZD) + (PD * CNZ) 655 : FRET = (UND * RHTD) 656 : C 657 : COEF = 0.5D0 * SURF / VOLG 658 : C 659 : BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG) 660 : & + FR * COEF 661 : BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG) 662 : & + FRUX * COEF 663 : BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG) 664 : & + FRUY * COEF 665 : BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG) 666 : & + FRUZ * COEF 667 : BN0.VAL(5,NLCEG) = BN0.VAL(5,NLCEG) 668 : & + FRET * COEF 669 : C 670 : C************* Viscous contribution 671 : C 672 : COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLG 673 : C 674 : BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG) 675 : & - (COEF * ROD) 676 : BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG) 677 : & - (COEF * RUXD) 678 : BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG) 679 : & - (COEF * RUYD) 680 : BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG) 681 : & - (COEF * RUZD) 682 : BN0.VAL(5,NLCEG) = BN0.VAL(5,NLCEG) 683 : & - (COEF * RETD) 684 : C 685 : C************* On calcule CN0 686 : C 687 : C CN0 in the referred report is given by 688 : C \sum_j 689 : C \left\{ 690 : C \frac{S_j}{2 \Omega_i} ~ 691 : C \left( 692 : C \sigma_{i,j}^{0} 693 : C \mathbf{U}_{j,R}^{n+1,m} 694 : C \right) 695 : C \right\} 696 : C 697 : C 698 : COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG) 699 : C 700 : CN0.VAL(1,NLCEG) = CN0.VAL(1,NLCEG) 701 : & + (COEF * ROD) 702 : CN0.VAL(2,NLCEG) = CN0.VAL(2,NLCEG) 703 : & + (COEF * RUXD) 704 : CN0.VAL(3,NLCEG) = CN0.VAL(3,NLCEG) 705 : & + (COEF * RUYD) 706 : CN0.VAL(4,NLCEG) = CN0.VAL(4,NLCEG) 707 : & + (COEF * RUZD) 708 : CN0.VAL(5,NLCEG) = CN0.VAL(5,NLCEG) 709 : & + (COEF * RETD) 710 : 711 : ENDDO 712 : C 713 : C********** End of loop on the faces of the element 714 : C 715 : C N.B.: ICEN = NLCEG 716 : C 717 : C 718 : C Computation of the low-Mach preconditioning Matrix 719 : C 720 : CG = (GAMG * PG / ROG)**0.5D0 721 : UCOG = MPOCO1.VPOCHA(ICEN,1) 722 : QG = 2.0D0 * RECG / 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 : Q2.VAL(1,ICEN) = (0.5D0 * QG * QG) 738 : Q2.VAL(2,ICEN) = -1.0D0 * RUXG / ROG 739 : Q2.VAL(3,ICEN) = -1.0D0 * RUYG / ROG 740 : Q2.VAL(4,ICEN) = -1.0D0 * RUZG / ROG 741 : Q2.VAL(5,ICEN) = 1.0D0 742 : C 743 : COEF1 = (GAMG - 1.0D0) / (CG * CG) 744 : COEF = (1.0D0 / COEF1) + (0.5D0 * QG * QG) 745 : C COEF=HT 746 : Q1.VAL(1,ICEN) = COEF1 747 : Q1.VAL(2,ICEN) = COEF1 * RUXG / ROG 748 : Q1.VAL(3,ICEN) = COEF1 * RUYG / ROG 749 : Q1.VAL(4,ICEN) = COEF1 * RUZG / ROG 750 : Q1.VAL(5,ICEN) = COEF1 * COEF 751 : C 752 : C 753 : C Computation of the diagonal coefficients 754 : C 755 : AN0.VALVEC(ICEN) = AN0.VALVEC(ICEN) 756 : & + (USDTI.VALVEC(ICEN) * 2.0D0 / DCFL) 757 : C 758 : DN0.VALVEC(ICEN) = DN0.VALVEC(ICEN) + (1.0D0 / DTPS) 759 : c & + (USDTI.VALVEC(ICEN) * 2.0D0 / DCFL) 760 : C 761 : ENDDO 762 : ENDDO 763 : C 764 : C*********************************************************************** 765 : C**** END FIRST JACOBI ITERATION *************************************** 766 : C*********************************************************************** 767 : C 768 : C*********************************************************************** 769 : C**** JACOBI LOOP ****************************************************** 770 : C*********************************************************************** 771 : C 772 : C 773 : LOGJAC=.TRUE. 774 : C 775 : DO IBLJAC=1,NJAC,1 776 : C 777 : C 778 : C Si ICOEF=2 (i.e. si on a "LJACOF" (voir kfm1.eso)) 779 : C alors relaxation SGS avec uniquement des balayages "forward" 780 : C 781 : C Si ICOEF=3 (i.e. si on a "LJACOB") 782 : C alors relaxation SGS avec uniquement des balayages "backward" 783 : C 784 : C Sinon (i.e. si on a "LJACOFB") 785 : C alors relaxation SGS avec des balayages "forward" et "backward" 786 : C selon la parité de la sous-itération en cours (IBLJAC)... 787 : C 788 : IF(ICOEF .EQ. 2)THEN 789 : LOGJAC=.TRUE. 790 : ELSEIF(ICOEF .EQ. 3)THEN 791 : LOGJAC=.FALSE. 792 : ELSE 793 : LOGINV=LOGJAC 794 : LOGJAC=(IBLJAC / 2 * 2) .EQ. IBLJAC 795 : IF(LOGJAC)THEN 796 : LOGJAC= .NOT. LOGINV 797 : ELSE 798 : LOGJAC= LOGINV 799 : ENDIF 800 : ENDIF 801 : C 802 : IF(LOGJAC)THEN 803 : ICEN=0 804 : IFACG=0 805 : ID=1 806 : ISOUP1=1 807 : ISOUP2=NSOUPO 808 : ELSE 809 : ICEN=NCEN+1 810 : IFACG=IFACT+1 811 : ID=-1 812 : ISOUP1=NSOUPO 813 : ISOUP2=1 814 : ENDIF 815 : C 816 : DO ISOUP=ISOUP1,ISOUP2,ID 817 : IPT2=MLESUP.LECT(ISOUP) 818 : NFACE=IPT2.NUM(/1) 819 : NELT=IPT2.NUM(/2) 820 : C 821 : IF(LOGJAC)THEN 822 : IELT1=1 823 : IELT2=NELT 824 : IFAC1=1 825 : IFAC2=NFACE 826 : ELSE 827 : IELT2=1 828 : IELT1=NELT 829 : IFAC2=1 830 : IFAC1=NFACE 831 : ENDIF 832 : C 833 : DO IELT=IELT1,IELT2,ID 834 : ICEN=ICEN+ID 835 : DO IFAC=IFAC1,IFAC2,ID 836 : C IFACG increase with IFAC,IELT,ISOUP 837 : IFACG = IFACG+ID 838 : NGCF = IPT2.NUM(IFAC,IELT) 839 : NLCF = MLENT2.LECT(NGCF) 840 : NGCEG = MELEFL.NUM(1,NLCF) 841 : NGCED = MELEFL.NUM(3,NLCF) 842 : NLCEG = MLENT1.LECT(NGCEG) 843 : NLCED = MLENT1.LECT(NGCED) 844 : C 845 : CNX = MPNORM.VPOCHA(NLCF,7) 846 : CNY = MPNORM.VPOCHA(NLCF,8) 847 : CNZ = MPNORM.VPOCHA(NLCF,9) 848 : CTX = MPNORM.VPOCHA(NLCF,1) 849 : CTY = MPNORM.VPOCHA(NLCF,2) 850 : CTZ = MPNORM.VPOCHA(NLCF,3) 851 : CVX = MPNORM.VPOCHA(NLCF,4) 852 : CVY = MPNORM.VPOCHA(NLCF,5) 853 : CVZ = MPNORM.VPOCHA(NLCF,6) 854 : C 855 : C**************** Left or right? 856 : C 857 : IF(NLCEG .EQ. ICEN)THEN 858 : LOGINV=.FALSE. 859 : ELSEIF(NLCED .EQ. ICEN)THEN 860 : LOGINV=.TRUE. 861 : ELSE 862 : WRITE(IOIMP,*) 'Anomalie dans kfm14.eso' 863 : WRITE(IOIMP,*) 'Problème de SPG' 864 : C 21 2 865 : C Données incompatibles 865 : CALL ERREUR(21) 866 : GOTO 9999 867 : ENDIF 868 : C**************** If right, inversion 869 : IF(LOGINV) THEN 870 : NLCED=NLCEG 871 : NLCEG=ICEN 872 : CNX = -1.0D0*CNX 873 : CNY = -1.0D0*CNY 874 : CNZ = -1.0D0*CNZ 875 : CTX = -1.0D0*CTX 876 : CTY = -1.0D0*CTY 877 : CTZ = -1.0D0*CTZ 878 : CVX = -1.0D0*CVX 879 : CVY = -1.0D0*CVY 880 : CVZ = -1.0D0*CVZ 881 : ENDIF 882 : C 883 : SURF = MPOVSU.VPOCHA(NLCF,1) 884 : SIGMAF = SIGMA0.VALVEC(NLCF) 885 : C 886 : C**************** Recuperation des Etats "gauche" et "droite" 887 : C 888 : ROG = CONS1.VAL(1,NLCEG) 889 : RUXG = CONS1.VAL(2,NLCEG) 890 : RUYG = CONS1.VAL(3,NLCEG) 891 : RUZG = CONS1.VAL(4,NLCEG) 892 : RETG = CONS1.VAL(5,NLCEG) 893 : GAMG = CONS1.VAL(6,NLCEG) 894 : C 895 : UXG = RUXG / ROG 896 : UYG = RUYG / ROG 897 : UZG = RUZG / ROG 898 : C 899 : UNG = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ) 900 : UTG = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ) 901 : UVG = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ) 902 : C 903 : RECG = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2) 904 : RECG = RECG / ROG 905 : PG = (GAMG - 1.0D0) * (RETG - RECG) 906 : C 907 : VOLG = MPOVOL.VPOCHA(NLCEG,1) 908 : C 909 : ROD = CONS1.VAL(1,NLCED) 910 : RUXD = CONS1.VAL(2,NLCED) 911 : RUYD = CONS1.VAL(3,NLCED) 912 : RUZD = CONS1.VAL(4,NLCED) 913 : RETD = CONS1.VAL(5,NLCED) 914 : GAMD = CONS1.VAL(6,NLCED) 915 : C 916 : UXD = RUXD / ROD 917 : UYD = RUYD / ROD 918 : UZD = RUZD / ROD 919 : C 920 : UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ) 921 : UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ) 922 : UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ) 923 : C 924 : RECD = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2) 925 : RECD = RECD / ROD 926 : PD = (GAMD - 1.0D0) * (RETD - RECD) 927 : C 928 : IF(NLCEG .EQ. NLCED)THEN 929 : C Murs ou condition limite 930 : NLLIM=MLELIM.LECT(NGCF) 931 : IF(NLLIM .EQ.0)THEN 932 : C Mur 933 : UND = -1.0D0 * UNG 934 : UTD = UTG 935 : UVD = UVG 936 : UXD = UND * CNX + UTD * CTX + UVD * CVX 937 : UYD = UND * CNY + UTD * CTY + UVD * CVY 938 : UZD = UND * CNZ + UTD * CTZ + UVD * CVZ 939 : C 940 : RUXD = ROD*UXD 941 : RUYD = ROD*UYD 942 : RUZD = ROD * UZD 943 : C 944 : ELSE 945 : ROD = MPLIM.VPOCHA(NLLIM,1) 946 : UXD = MPLIM.VPOCHA(NLLIM,2) 947 : UYD = MPLIM.VPOCHA(NLLIM,3) 948 : UZD = MPLIM.VPOCHA(NLLIM,4) 949 : PD = MPLIM.VPOCHA(NLLIM,5) 950 : C 951 : RUXD = ROD*UXD 952 : RUYD = ROD*UYD 953 : RUZD = ROD * UZD 954 : C 955 : GAMD = GAMG 956 : C 957 : UND = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ) 958 : UTD = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ) 959 : UVD = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ) 960 : RETD = ((1.0D0/(GAMD - 1.0D0))*PD) 961 : & +(0.5D0*ROD*(UXD**2+UYD**2+UZD**2)) 962 : ENDIF 963 : ENDIF 964 : C 965 : C**************** We check that density and pressure are positive 966 : C 967 : IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR. 968 : & (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN 969 : WRITE(IOIMP,*) 'FMM' 970 : WRITE(IOIMP,*) 'ITERATION JAC',IBLJAC 971 : WRITE(IOIMP,*) 'ITERATION ICEN',ICEN 972 : WRITE(IOIMP,*) ROG,PG 973 : WRITE(IOIMP,*) ROD,PD 974 : C 975 : C 22 0 976 : C******************* Opération malvenue. Résultat douteux 977 : C 977 : CALL ERREUR(22) 978 : IPROBL = 1 979 : GOTO 9998 980 : ENDIF 981 : C 982 : C 983 : C**************** We compute the flux 984 : C FLU.VALVEC(1:4,IFAC) flux in IFAC face 985 : C 1 -> mass 986 : C 2 -> momentum (x) 987 : C 3 -> momentum (y) 988 : C 4 -> energy 989 : C 990 : RHTD = RETD + PD 991 : FR = (ROD * UND) 992 : FRUX = (ROD * UND * UXD) + (PD * CNX) 993 : FRUY = (ROD * UND * UYD) + (PD * CNY) 994 : FRUZ = (ROD * UND * UZD) + (PD * CNZ) 995 : FRET = (UND * RHTD) 996 : C 997 : COEF = 0.5D0 * SURF / VOLG 998 : C 999 : BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) 1000 : & + FR * COEF 1001 : BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) 1002 : & + FRUX * COEF 1003 : BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) 1004 : & + FRUY * COEF 1005 : BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) 1006 : & + FRUZ * COEF 1007 : BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG) 1008 : & + FRET * COEF 1009 : C 1010 : C**************** Viscous contribution 1011 : C 1012 : COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLG 1013 : C 1014 : BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) 1015 : & - (COEF * ROD) 1016 : BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) 1017 : & - (COEF * RUXD) 1018 : BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) 1019 : & - (COEF * RUYD) 1020 : BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) 1021 : & - (COEF * RUZD) 1022 : BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG) 1023 : & - (COEF * RETD) 1024 : C 1025 : C**************** On calcule DIS1 1026 : C 1027 : C DIS1 in the referred report is given by 1028 : C \sum_j 1029 : C \left\{ 1030 : C \frac{S_j}{2 \Omega_i} ~ 1031 : C \left( 1032 : C \sigma_{i,j}^{0} 1033 : C \mathbf{U}_{j,R}^{n+1,m,l} 1034 : C \right) 1035 : C \right\} 1036 : C 1037 : C 1038 : COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG) 1039 : C 1040 : CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG) 1041 : & + (COEF * ROD) 1042 : CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG) 1043 : & + (COEF * RUXD) 1044 : CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG) 1045 : & + (COEF * RUYD) 1046 : CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG) 1047 : & + (COEF * RUZD) 1048 : CN1.VAL(5,NLCEG) = CN1.VAL(5,NLCEG) 1049 : & + (COEF * RETD) 1050 : 1051 : ENDDO 1052 : C 1053 : C************* End of the loop on the faces of the element 1054 : C 1055 : DO ICOMP=1,5,1 1056 : D1.VALVEC(ICOMP) = 0.0D0 1057 : ENDDO 1058 : C 1059 : COEF=0.0D0 1060 : DO I1=1,5,1 1061 : COEF = COEF 1062 : & + (Q2.VAL(I1,ICEN) * (CN1.VAL(I1,ICEN) 1063 : & -CN0.VAL(I1,ICEN))) 1064 : ENDDO 1065 : COEF = COEF 1066 : & * (1.0D0 - PHI2V.VALVEC(ICEN)) 1067 : & / PHI2V.VALVEC(ICEN) 1068 : C 1069 : DO I1=1,5,1 1070 : D1.VALVEC(I1) = (CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN)) 1071 : & + (COEF * Q1.VAL(I1,ICEN)) 1072 : ENDDO 1073 : C 1074 : C************* At this moment D1 = \Gamma \cdot (\Delta CN1) 1075 : C 1076 : DO ICOMP=1,5,1 1077 : C 1078 : C************* First increment 1079 : C 1080 : D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP) 1081 : & + RES.VAL(ICOMP,ICEN) 1082 : & + (BN0.VAL(ICOMP,ICEN) 1083 : & -BN1.VAL(ICOMP,ICEN))) 1084 : & / (AN0.VALVEC(ICEN) 1085 : & +DN0.VALVEC(ICEN)) 1086 : ENDDO 1087 : C 1088 : C************* At this moment, in the referred report, 1089 : C D1 is (\delta \mathbf{U}')^{l+1}_i in section A3 1090 : C 1091 : C************* Second increment 1092 : C 1093 : COEF=0.0D0 1094 : DO I1=1,5,1 1095 : COEF = COEF + (Q2.VAL(I1,ICEN) * D1.VALVEC(I1)) 1096 : ENDDO 1097 : C 1098 : C************* COEF is the scalar product between Q2 and the first 1099 : C increment. 1100 : C Let us multiply it by 1101 : C b_i^0 / (a_i^0 + b_i^0) 1102 : C 1103 : COEF = COEF 1104 : & * AN0.VALVEC(ICEN) 1105 : & * (PHI2V.VALVEC(ICEN) - 1.0D0) 1106 : COEF = COEF 1107 : & / (AN0.VALVEC(ICEN) 1108 : & + (DN0.VALVEC(ICEN) * PHI2V.VALVEC(ICEN))) 1109 : C 1110 : DO I1=1,5,1 1111 : D1.VALVEC(I1) = D1.VALVEC(I1) 1112 : & + COEF * Q1.VAL(I1,ICEN) 1113 : ENDDO 1114 : C 1115 : C************* D1 = \delta \mathbf{U} 1116 : C 1117 : DO ICOMP=1,5,1 1118 : DU.VAL(ICOMP,ICEN) = D1.VALVEC(ICOMP) 1119 : C 1120 : CONS1.VAL(ICOMP,ICEN) = CONS0.VAL(ICOMP,ICEN) 1121 : & + D1.VALVEC(ICOMP) 1122 : C 1123 : BN1.VAL(ICOMP,ICEN) = 0.0D0 1124 : CN1.VAL(ICOMP,ICEN) = 0.0D0 1125 : ENDDO 1126 : C 1127 : C********** End of the loop on the elements 1128 : ENDDO 1129 : C 1130 : C******* End of the loop on the Domains 1131 : ENDDO 1132 : C 1133 : C**** End of the Jacobi loop 1134 : ENDDO 1135 : C 1136 : 9998 CONTINUE 1137 : DO ICEN=1,NCEN,1 1138 : DO ICOMP=1,5,1 1139 : MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN) 1140 : ENDDO 1141 : ENDDO 1142 : C 1143 : SEGSUP DU 1144 : SEGSUP RES 1145 : SEGSUP CONS0 1146 : SEGSUP CONS1 1147 : SEGSUP Q1 1148 : SEGSUP Q2 1149 : SEGSUP SIGMA0 1150 : SEGSUP AN0 1151 : SEGSUP AN0 1152 : SEGSUP BN0 1153 : SEGSUP BN1 1154 : SEGSUP CN0 1155 : SEGSUP CN1 1156 : SEGSUP DN0 1157 : SEGSUP PHI2V 1158 : SEGSUP USDTI 1159 : C 1160 : SEGDES MPOVSU 1161 : SEGDES MPOVDI 1162 : SEGDES MPOVOL 1163 : SEGDES MPNORM 1164 : SEGDES MPOCO 1165 : SEGDES MPCOEV 1166 : C 1167 : SEGSUP MPOCO1 1168 : C 1169 : SEGDES MPODU 1170 : C 1171 : SEGDES MELEFL 1172 : SEGSUP MLENT1 1173 : SEGSUP MLENT2 1174 : C 1175 : SEGSUP MLELIM 1176 : IF(MPLIM .GT. 0) SEGDES MPLIM 1177 : C 1178 : DO ISOUP=1,NSOUPO,1 1179 : IPT2=MLESUP.LECT(ISOUP) 1180 : SEGDES IPT2 1181 : ENDDO 1182 : IPT2 = MELTFA 1183 : IF(ISOUP .GT. 1) SEGDES IPT2 1184 : SEGSUP MLESUP 1185 : C 1186 : 9999 CONTINUE 1187 : RETURN 1188 : END 1189 : C 1190 : 1191 : 1192 :
© Cast3M 2003 - All rights reserved.
Disclaimer